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Summary 


A computational technique is developed that is suitable for performing prelimi- 
nary design aeroelastic and structural dynamic analyses of large aspect ratio lifting 
surfaces. The method proves to be quite general and can be adapted to solving 
various t^vo-point boundary value problems. 

The solution method, which is applicable to both fixed and rotating wing 
configurations, is based upon a formulation of the structural equilibrium equations 
in terms of a hybrid state vector containing generalized force and displacement 
variables. A mixed variational formulation is presented that conveniently yields a 
useful form for these state vector differential equations. Solutions to these equations 
are obtained by employing an integrating matrix method. The application of an 
integrating matrix provides a discretization of the differential equations that only 
requires solutions of standard linear matrix systems. It is demonstrated that matrix 
partitioning can be used to reduce the order of the required solutions. Results are 
presented for several example problems in structural dynamics and aeroelasticity to 
verify the technique and to demonstrate its use. These problems examine various 
types of loading and boundary conditions and include aeroelastic analyses of lifting 
surfaces constructed from anisotropic composite materials. 

Integrating matrices, which provide a powerful tool for solving differential 
equations, are discussed in detail, and methods are given for their calculation. 
A derivation and calculation procedure is presented for a new type of maximum 
accuracy integrating matrix based upon orthogonal polynomials. 
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Chapter 1 


Introduction 


An expanded utilization of laminated composite materials in aircraft 
structural design has led to a search for new ways to employ these relatively high 
specific strength and stiffness materials. One result of this search has been the 
development of the concept of aeroelastic tailoring of a lifting surface, in which the 
directional characteristics of the composite material are used to synthesize a struc- 
ture with enhanced aeroelastic performance. But along with the possibility for in- 
novative design with structural composites comes greater complexity in the prelimi- 
nary design task. This increased complexity arises in part from the anisotropic 
nature of the composites materials and in part from the increased design freedom 
allowed by these materials. Because of the additional complexity of the design task, 
new analysis tools are needed to aid the preliminary designer in efficiently evaluat- 
ing the sometimes large number of design concepts available to him. Therefore, 
the primary objective of this research has been to develop a simple and versatile 
analysis method compatible with the needs of preliminary aeroelastic and structural 
dynamic design. 

The motivation for this research effort stemmed from a desire to investigate the 
performance enhancements that can be achieved by aeroelastically tailoring large 
aspect ratio composite lifting surfaces. Essentially, aeroelastic tailoring involves 
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designing a structure to take advantage of the elastic deformation during loading. 
For static aeroelastic problems, this means controlling the relative amounts of bend- 
ing and torsional deflection of a wing or lifting surface. By maintaining a desirable 
wing deformation shape, or by passively controlling the distribution of aerodynamic 
loading, it is often possible to enhance aerodynamic performance and to extend 
the operating envelopes of a lifting surface. For dynamic aeroelastic problems, the 
coupling between bending and torsion of composite structures provides a way of 
maximizing the dynamic instability (flutter) speed of a lifting surface. Since the 
primary objective of this research is to develop a convenient method for analyzing 
such aeroelastic and dynamic behavior, the above mentioned problems provide some 
excellent, nontrivial examples for verification of the devised solution method. At the 
same time, these example solutions hopefully provide a firm foundation for other 
in-depth studies of the aeroelastic behavior of composite structures, including the 
investigation of optimized aeroelastic designs. 

Historically, most aeroelastic analyses of composite structures have been carried 
out by very complex computer codes involving finite element structural methods 
coupled with lifting surface aerodynamics. Unfortunately, these complicated numer- 
ical approaches can tend to obscure a basic understanding of the important param- 
eters appearing in the analysis and, owing to cost considerations, often preclude 
an extensive study involving numerous design variations. Recent developments, 
such as those of Gimmestad [1], offer a suitable alternative for preliminary design 
investigations. 

A fundamental approach to performing the aeroelastic and dynamic analyses of 
a structure described by one independent spatial coordinate involves formulating the 
ordinary differential equations representing the aeroelastic or dynamic response and 
obtaining analytical solutions to the resulting boundary value problems. Although 
the coupled bending and torsion equations can be formulated, it is often difficult, 
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if not impossible, to obtain analytical solutions for the general case in which the 
coefficients of the equations are variable. Some useful solutions have been obtained, 
however, for cases in which the coefficients in the linear aeroelastic equations can 
be written as constants. For instance, solutions to the differential equations for 
divergence and load distribution have been obtained for isotropic metallic wings 
by Diederich and Budiansky [2] and Diederich and Foss [3], and more recently, 
divergence and load distribution solutions for composite swept forward wings have 
been obtained by Weisshaar [4,5]. But with the application of the hybrid state vector 
approach discussed herein, approximate solutions to the differential equations can 
be easily obtained for much more difficult cases involving variable coefficients. The 
hybrid state vector approach has been utilized by Lehman [6] to obtain a variety 
of aeroelastic solutions, including flutter of composite wings. This type of solution 
does not require an explicit calculation of structural influence coefficients and can 
utilize various forms of aerodynamic influence matrices. 

A major requirement for a solution method to be used in preliminary design 
is that the method be reasonably flexible in allowing solution of different types 
of problems, and yet easily specialized so that computations can be carried out 
efficiently. Furthermore, it is desirable to have a numerical solution that is easily 
programmable and that makes use of standard numerical methods, thus requiring 
minimal investment in software. It has been found that these requirements are 
well satisfied by a mixed state vector formulation of the differential equations 
combined with an integrating matrix solution procedure — hence, one of the reasons 
for referring to the method as a hybrid approach. 

Other investigators (see the introduction to Chapter 3) have provided initial for- 
mulations for the concept of the integrating matrix solution in structural mechanics 
and have applied this concept to solving a variety of problems. Compared to other 
numerical approaches, such as finite element and finite difference, relatively little 
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has been done to generalize the integrating matrix method. In order to bring some 
generality to the integrating matrix method, it has been found useful to incorporate 
some of the familiar concepts employed in finite element analysis. 

Even though the solution method presented here is applied only to two-point 
boundary value problems arising in aeroelasticity and structural dynamics, the 
approach is, in fact, quite general and can be applied to initial value problems as 
well as to multipoint boundary value problems. The method can also be extended 
to handle systems described by more than one independent variable. 

The compact matrix notation used in the development of the hybrid state vector 
method is intended to aid in the task of programming the solutions, regardless of 
the programming language employed. Solutions for dynamic flutter instabilities 
in Chapter 7, which are iterative by nature, operate quite efficiently in languages 
like FORTRAN or Pascal. The hybrid state vector solution formulations, however, 
are especially suited to the matrix oriented programming language APL. In fact, 
the hybrid state vector method presented here, when coupled with APL, forms an 
extremely powerful interactive problem solving tool. It is further anticipated that 
the hybrid state vector solutions, since they are formulated in terms of simultaneous 
matrix operations, will be readily adaptable to parallel processing techniques. 

In Chapter 2, a mixed variational formulation is presented for obtaining the 
linear state vector differential equations of structural equilibrium. This formulation 
is given for structures that can be described by one independent spatial variable. 
By casting the aerodynamic and inertial loads acting on a structure in terms of 
the displacement state variables, the state vector equations can be expanded into 
the equations suitable for aeroelastic and structural dynamic analyses. A detailed 
form of the state vector equations is then presented to describe an anisotropic plate- 
beam that is constructed from laminated composite materials. These equations are 
reserved for later use in the example solutions of Chapters 6 and 7. 
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Chapter 3 gives a general derivation of integrating matrices and describes how 
they are applied to the integration of either continuous or discontinuous integrands. 
Several different types of integrating matrices are discussed, including maximum 
precision integrating matrices based upon orthogonal polynomial approximations. 
The concept of a differentiating matrix is also introduced. 

Chapter 4 describes how integrating matrices are used to formulate solutions 
for the discretized versions of the state vector equations derived in Chapter 2. By 
using matrix partitioning techniques, it is then shown that reduced order matrix 
equations for the displacement variables can be obtained by eliminating the force 
variables. 

In Chapter 5, sample solutions are presented for simple beam and rod problems. 
These examples illustrate the application of the hybrid state vector method to 
the solution of two-point boundary value problems. Continuous and discontinuous 
parameter problems are demonstrated along with various types of boundary and 
loading conditions. Numerical results are compared with analytical results to 
evaluate the accuracy of the integrating matrix solutions. 

Chapter 6 presents sample solutions for divergence and elastic lift distribution 
of composite wings. For the composite wings, solutions are given for the case of 
forward aerodynamic sweep. Brief comparisons are made with alternate solutions 
available for these problems. 

Chapter 7 demonstrates solutions for flutter instabilities of both isotropic and 
composite wings. The isotropic wing flutter solutions are compared with known 
analytical solutions. 

A brief summary and recommendations for additional research are given in 
Chapter 8. 
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Chapter 2 


Hybrid State Vector Equations 


A GENERAL FORMULATION of the structural equations will be presented 
which casts them into a state vector form involving a coupled system of first order 
differential equations. The state vector that appears in the following derivations 
will be termed a hybrid state vector in the sense that it is derived from a so-called 
mixed formulation involving both stress and displacement variables. Although this 
formulation is not entirely new to structural mechanics, it has not seen extensive 
use, nor has it been included among the everyday tools of most engineers work- 
ing in structures and structural dynamics. However, some very noteworthy de- 
velopments of improved numerical procedures based on mixed formulations com- 
bined with finite-difference solutions have been reported by Noor, Stephens, and 
Fulton [7]. Additional work presented by Noor and Stephens [8,9] has further 
demonstrated both the simplicity and high accuracy of such mixed formulation 
procedures. Results obtained by Stroud and Mayers [10] indicate that a numerical 
solution based upon direct application of a mixed variational principle also offers 
superior accuracy and convergence, especially for bending-moment solutions. More 
recently, investigations by Steele [11], Steele et al. [12], and Steele and Barry [13] 
have indicated that mixed state vector formulations of the differential equations 
in conjunction with asymptotic solutions can be advantageous for both analytical 
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investigation and numerical calculation. The present investigation will demonstrate 
that a simple and highly accurate numerical solution procedure is obtained by com- 
bining the mixed state vector formulation of the structural differential equations 
with an integrating matrix solution approach. It is worth noting that the transfer 
matrix method of structural solution also employs a mixed state vector similar to 
the one used in the following formulations. 

The derivations to be presented employ a mixed variational formulation that 
can be consistently applied to a broad class of structural problems. This formula- 
tion, which will be discussed in Section 2.1, involves terms that are expressible as 
a product of generalized stresses and strains in addition to other terms that can be 
related to the complementary energy. For systems with linear stress-strain behavior, 
the complementary energy is, of course, the same as the strain energy. Therefore, 
in the context of linear systems, this formulation is equivalent to the more usual 
stationary potential energy approach. As will be demonstrated, the mixed varia- 
tional formulation provides a convenient way of expressing the energy functional 
and allows a direct determination of the state vector equations in a desirable form. 

There are also many classical structural problems for which differential equa- 
tions already exist. In these instances, it may be convenient to recast these equations 
into a matrix form directly and dispense with the formality of rederiving them. As 
many readers are well aware, it is possible to take higher order differential equations 
and convert them to an equivalent system of first order equations. But this process 
becomes increasingly difficult as the complexity of the system increases. Regardless 
of how one chooses to obtain the differential equations describing a structural prob- 
lem, there is a preferred way to write the state vector equivalent. The preferred 
state vector form will be shown to arise naturally from a mixed variational for- 
mulation. As will be discussed in more detail in later sections, the mixed (hybrid) 
state vector form of the equations, with fundamental unknowns consisting of both 
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generalized force and displacement parameters, will possesses a natural decomposi- 
tion that proves quite useful for numerical solutions. 

The equation derivations that follow will be presented in matrix notation. The 
primary advantage of matrix notation is that it allows a consistent treatment of 
problems with differing size and complexity. Furthermore, it is anticipated that 
the equation formulation will lead directly to a numerical algorithm that must 
necessarily deal with the matrix terminology of discretized systems. An additional 
reason for adhering to matrix notation is that several intermediate points exist at 
which further analytical formulation and simplification can be put aside in favor of 
numerical computation. If a matrix formulation is used throughout, it is easy to 
begin numerical calculations at these intermediate points. 

2.1 Variational Formulation of State Vector Equations 

A mixed variational formulation is presented here for the hybrid state vector 
equations that describe structural equilibrium. For the interested reader, some 
rather general examples of mixed variational statements can be found, for instance, 
in Nemat-Nasser [14], Also, brief historical accounts of mixed variational methods 
in solid mechanics appear in both Nemat-Nasser [15] and Reissner [16]. In these 
accounts, the work of a number of investigators is cited, including the work of 
authors such as Hellinger, Reissner, and Washizu. In the literature in general, the 
mixed variational formulation involving independent variation of both stress and 
displacement variables is usually referred to as a Reissner (or sometimes Hellinger- 
Reissner) formulation, whereas the principle involving independent variation of 
stress, displacement, and strain is often referred to as a Hellinger-Reissner-Washizu 
formulation. In the work presented here, the development of the mixed state 
vector equations is based upon the formulation of Reissner [17] in which stress and 
displacement variables are independently varied to yield the appropriate equations 


- 11 - 



and boundary conditions. A similar development of mixed state vector equations 
from the Reissner formulation is given in Ref. [11]. 

Although the Reissner formulation is applicable to either linear or nonlinear 
problems, the following presentation will restrict consideration to the solution of 
linear aeroelastic and structural dynamic equations. From nonlinear equations, it 
is often feasible to obtain a set of linearized equations by perturbing about an 
appropriate nonlinear solution. This might be useful, for example, when consider- 
ing problems with geometric nonlinearities. The linear perturbation equations ob- 
tained by such an approach fall within the scope of the following linear analyses. 
When considering future extensions of the present work to nonlinear analyses, it 
is anticipated that the Reissner principle should prove valuable for problems with 
nonlinear material behavior. Some important illustrations of the application of the 
Reissner principle to problems involving nonlinear material behavior are given by 
Nimmer and Mayers [18] and Anderson and Mayers [19]. 

For those problems that can be described by a single spatial variable, x, the 
mixed variational formulation can be written in general terms as 


8 



( 2 . 1 ) 


where the prime on y indicates partial differentiation with respect to x only. The 
Euler-Lagrange equations resulting from variation on x are 


dx dy' 


du 

dy 


( 2 . 2 ) 


Since the variation is being taken only with respect to the spatial variable, the time 
variable simply follows along as a parameter. For static problems, time disappears 
from the previous equations. 

The next step in the formulation is to give an appropriate form of the functional 
appearing in Eqs. (2. 1-2.2). When linearity is invoked, it is then possible to express 
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the functional in the following convenient form: 


i/=y'^Dy-iy^Ky + p^y (2.3) 

where K is a spatially variable symmetric matrix containing structural relations, 
p contains the external loads, and D is defined such that (D — D^) is a constant 
skew-symmetric matrix with unit elements. 

The state vector, y, is specified in the form 

y == {y/- VdV (2-4) 

where y/^ are generalized forces and y/> are generalized displacements. This form 
for y is the same as would appear in a transfer matrix structural solution based 
upon a “mixed” finite element force-displacement relationship. The precise nature 
of the matrix terms appearing in the above representation of the functional will be 
clarified with specific examples. 

Next, substituting Eq. (2.3) into Eq. (2.2) and performing the indicated differen- 
tiation yields the equilibrium equations 

-Jy'-Ky + p = 0 (2.5) 


in which 


Noticing that 


J = D - = 



J ^ = — J and hence J^J == I , 


Eq. (2.5) can be rewritten in the standard state vector form 


y' = Zy - a 


( 2 . 6 ) 


(2.7) 


( 2 . 8 ) 
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where 


Z = JK and a = Jp. (2.9) 

Eq. (2.8) will be the starting point for analyses to be presented in later chapters. 

Rather than using the procedure demonstrated above, Eq. (2.5) and correspond- 
ing boundary conditions can also be obtained by substituting Eq. (2.3) into Eq. (2.1) 
and applying integration by parts. The consistent boundary conditions obtained 
with this approach are one of the key advantages of a variational development. 
Corresponding to the state vector equations in Eq. (2.8), the boundary conditions 
are obtained as 

y^D^(^y)lo = 0- (210) 

Fortunately, in the state vector formulation, these boundary conditions always 
remain quite simple. This will prove to be especially advantageous when dealing 
with anisotropic structures, for which other formulations can yield coupled and 
considerably more complicated forms of the boundary conditions. 

At this point, a simple example describing the lateral bending deflection of a 
beam will help to clarify the nature of the matrix terms appearing in the foregoing 
derivation. The example presented here follows an example given by Steele [11] for 
a Timoshenko beam. In the notation used in this study, the functional in Eq. (2.3) 
can be written for a Timoshenko beam in the form 

U - - pw ( 2 . 11 ) 

where the moment resultant is Mx', the transverse shear resultant is Vx; the rotation 
of the normal is 7; the normal displacement is w; the shear compliance is ^e] and 
the load per unit length is p. Prime, of course, denotes differentiation with respect 
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to X. In order to obtain a resemblance to the expression in Eq. (2.3), Eq. (2.11) can 
be rewritten as 


u = - i((E/)-*M| + - pw- (2.12) 

For linear systems, this functional can be constructed by subtracting the com- 
plementary energy from the product of stresses and strains, after which a slight 
rearrangement yields Eq. (2.12). Also note that the nature of U, as used in this 
presentation, implies that it has been obtained by integrating an energy density 
functional over the cross section of the beam. 

If the state vector is now defined as 

y={Mr 7 (2.13) 

then p, D, and K in Eq. (2.3) take on the forms 

p={0 0 0 -p}"^ (2.14) 


D = 


■0 

0 

1 

.0 


0 0 0 ' 

0 0 0 

0 0 0 

10 0 . 


(2.15) 


0 0 O' 

0 ^,-10 

0 -10 0 

- 0 0 0 0 - 


(2.16) 


The procedure for arriving at this form of the matrix equations first involves 
specifying the state vector, y, which is taken to be the same vector as would be 
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obtained in a mixed, or hybrid force-displacement formulation of the equilibrium 
equations. Such a formulation is discussed, for example, in Chapter 2 of McGuire 
and Gallagher [20]. This state vector form is also the same as that used in the 
transfer matrix format of structural analysis. It is assumed here, as shown in Eq. 
(2.4), that y will always be partitioned into two sets; one set contains generalized 
forces and the other contains generalized displacements. Once the state vector is 
specified in this way, then p is chosen so that one obtains the proper potential of 
the external loads. 

After choosing the state vector, then D must be determined such that the 
derivative terms in the functional are given by the first group on the right hand 
side of Eq. (2.3). This requirement is met by specifying that D always be a square 
matrix having the same form as that given in Eq. (2.15); that is, it should have the 
same structure as Eq. (2.15) and should always contain only zero and unit terms. 
One can refer to Eq. (2.6) to see that D must be specified in this way to insure that 
J be an antisymmetric matrix with unit elements. If J is as shown in Eq. (2.6), 
then Eq. (2.5) and Eq. (2.8) are said to have a symplectic character. The symplectic 
nature of these equations means that an especially simple relationship will exist 
between the fundamental solution of the system in Eq. (2.8) and its adjoint (see 
Eq. (G.17) in Appendix G). For a description of this useful property of symplectic 
systems, refer to page 157 of Bryson and Ho [21]. 

Other remaining terms in the functional are now determined by specifying K, 
which is restricted to be a symmetric matrix. The elements of K contain spatially 
dependent constitutive terms and fixed structural kinematic relationships. It is 
usually easy to determine the elements of K by simple observation. 

Having determined y, p, D, and K as just discussed, it is then clear that a and 
Z in Eq. (2.9) become 
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a = {0 p 0 0}^, 


( 2 . 17 ) 


r 0 1 0 01 


z = 


0 


0 

0 


L 0 0 , 


0 0 

0 0 

-1 0 . 


( 2 . 18 ) 


It should be noted that the right hand side of Eq. (2.18) displays a particular 
form that will show up again, even for much more complicated problems than 
the one discussed above. As mentioned in the introduction to this chapter, for 
many simple problems it is not necessary to follow through the complete equation 
derivation as just presented. In fact, once one becomes familiar with the matrix 
formulation given here, it is usually easy to write down Eq. (2.18) directly, at 
least for relatively simple systems. However, for considerably more complicated 
situations, the variational approach provides a consistent method for formulating 
equilibrium equations and boundary conditions. In addition, the mixed variational 
formulation also demonstrates a “natural” form of the state vector equations. Later, 
Section 4.2 will show that this natural form, when coupled with an integrating 
matrix solution procedure, provides significant savings in the numerical solution 
by allowing convenient analytical simplification. One should also keep in mind 
that with knowledge of the natural form of the equations it is possible to recast 
equations derived by other methods. To give one example, the nonlinear equations 
of an initially bent and twisted rod derived in Chapter 18 of Love [22] and examined 
by Ojalvo and Newman [23] can be linearized and recast into the desired form. In 
fact, a similar approach has been followed by Nitzsche [24] to obtain hybrid state 
vector equations for the aeroelastic analysis of vertical axis wind turbines. 
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2.2 Aeroelastic and Structural Dynamic Equations 


The general equations presented in the previous section can be specialized to 
both aeroelastic and structural dynamic problems by being more specific about the 
type of loading. For both aeroelastic and dynamic problems the loading can be 
related in some way to the displacements of the structure (i.e., the displacement 
states of the structural state vector). In static aeroelastic problems the airloads are 
directly determined by the deformed shape of the structure, whereas in dynamic 
problems, the inertial, aerodynamic, and structural damping loads are related to 
time rates of change of the structural deformation. The discussion in this section 
will focus on the way in which these loads appear in the structural equilibrium 
equations, and the form of the equations to be used in later analyses will be given. 

First, static aeroelastic problems will be analyzed by the usual procedure of 
breaking the total external loads acting on the system into a summation of those 
loads that act on a rigid structure plus perturbation loads due to elastic deformation. 
Therefore, the loads vector a., which first appeared in Eq. (2.8), can be rewritten as 

a = gAy + 8r (2.19) 

where ar is a vector of the nonhomogeneous loads acting on a rigid structure, 
q is the dynamic pressure, and A is developed from an aerodynamic influence 
relationship. In fact, for the discrete version of these equations, A contains terms 
from an inverse aerodynamic influence matrix. Additionally, if A is partitioned 
corresponding to the force and displacement subsets of the state vector, only one of 
its submatrices contains nonzero elements, namely, that submatrix providing forces 
due to displacement. Although one can make use of various aerodynamic theories 
to calculate A, the analyses of Chapters 5-7 will primarily use aerodynamic strip 
theory. 
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The static aeroelastic equations are now obtained by substituting Eq. (2.19) 
into Eq. (2.8). The result is 

y' = Zy-gAy-ar, (2.20) 

with y = y(a:). When the dynamic pressure is specified, then Eq. (2.20) is simply a 
system of equations for y. On the other hand, if the nonhomogeneous term ar is 
set to zero and the dynamic pressure is left unspecified, then Eq. (2.20) leads to a 
divergence eigenvalue problem, with q being the divergence dynamic pressure. 

Now consider a dynamic aeroelastic system. The airloads still depend in some 
way on the displacements, but now time has entered the equations. Furthermore, 
inertia loads, and possibly damping loads, must be included in the analysis; as 
mentioned earlier these are also related in some way to the displacement. The 
approach taken here will be to remove the differential time dependence of the 
dynamic equations by Laplace transformation on time, thereby obtaining equations 
with an algebraic dependence on the Laplace variable, s. This has an added 
advantage for unsteady aeroelastic problems since the unsteady aerodynamic terms 
for general motion are conveniently described in the Laplace domain. 

After Laplace transformation, one finds that the homogeneous state vector 
equations suitable for aeroelastic stability analysis can be written in the form 

-^y = Zy + s^My + aCy — Q(a, g)y (2.21) 

ax 

where y = y(®, a) and the matrices M, C, and Q contain, respectively, the mass, 
damping, and unsteady aerodynamic terms. (The hat symbol denotes a Laplace 
transformed variable). As mentioned for the static aeroelastic problem in Eq. (2.20), 
the matrices expressing any form of displacement dependent loading (here, M, C, 
and Q) have only one nonzero partition, namely, the partition that multiplies the 
displacements in the state vector. 
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In order to specialize Eq. (2.21) to free vibration analysis, neglect the damping 
and aerodynamic terms, C and Q. For undamped vibration, the Laplace variable 
s is purely imaginary. Therefore, if one takes s = jcj, then the free vibration 
counterpart of Eq. (2.21) is given by 


—y = Zy - w^My. 


( 2 . 22 ) 


Clearly, dynamic response problems in the time domain can also be accommodated 
by including forcing terms on the right hand sides of the time domain versions of 
Eqs. (2.21) and (2.22). 


2.3 Anisotropic Beam Equations 

A simplified anisotropic plate-beam model is presented here for which the 
resulting equations are suitable for analyzing aeroelastic phenomena of large aspect 
ratio lifting surfaces. The purpose for developing these equations is twofold: first, 
they will help clarify the application of the foregoing general formulation and 
second, these equations will later be used in numerical examples of aeroelastic 
analyses. The assumptions used in developing the equations for the plate-beam 
model will be discussed briefly here, but it should be noted that they are similar to 
those employed by Weisshaar [4,5] to describe laminated composite box beam lifting 
surfaces. As a consequence, the following equations will be specialized somewhat to 
deal with structures whose anisotropic behavior arises due to laminated composite 
construction. A summary of composite plate lamination theory is presented in 
Appendix C; for a more detailed development, however, the reader can refer to 
Chapter 4 of the text by Jones [25]. 

With the aid of the cartesian coordinate system presented for the lifting surface 
model in Fig. 1, the aeroelastic equations will be developed for aerodynamic strip 
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sections taken normal to the structural reference axis. The structural reference 
axis is specified to be coincident with the x-axis, also shown in Fig. 1. Primarily for 
reasons of convenience, it will be assumed that this reference axis coincides with the 
geometric middle surface of the structural box. Although the reference axis location 
is arbitrary, this particular choice follows the conventions used for development of 
the composite plate constitutive relations. For this study, it is assumed that the 
structural reference axis is a straight tine which, of course, will be swept accordingly 
as the wing aerodynamic sweep angle, A, changes. When considering aerodynamic 
surfaces with structural axis curvature, then the present approach should be adapted 
to take this curvature into account. It is further assumed that no appreciable 
chordwise deformation occurs in cross sections normal to the structural reference 
axis, so that wing deformation is only a function of the spanwise coordinate, x. 
Applying these assumptions means that the deformation of the plate-beam model 
can be represented in terms of a bending deflection w(y), positive downward, of the 
reference axis, plus a rotation a(y), positive nose-up, about this axis. For problems 
dealing with rotating beams, an additional displacement variable u, along the x-axis, 
must be added to describe the deformation. In the presentation given here, bending 
deformation in the x-y plane is neglected. This lead-lag deformation can be easily 
added, however, in a more detailed analysis. 

As the next step, the in-plane strains and curvatures can be written in terms 
of the foregoing displacement variables by applying the differential relationships 
that describe strain-displacement for a plate (see Appendix C). But first, an addi- 
tional remark concerning shearing deformations should be made. Since standard 
lamination theory assumes no transverse shearing effects, this assumption will be 
adhered to here, but only for the composite laminate. That is, the shearing defor- 
mation of the laminated portion of the structure will be assumed negligible, but the 
gross shearing behavior of the overall structure can still be included. For instance. 
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Fig. 1. Lifting Surface Model 

standard composite structure fabrication quite often employs thin laminated face 
sheets placed over shear webs or a thick deformable core material. The transverse 
shearing eflfects induced by the webs or the core material can be included by adding 
a shearing energy to the energy functional. The transverse shear effects included 
here are assumed to arise in this manner. At the same time, it can reasonably be 
assumed that the composite cover sheets carry most of the bending stresses (and 
bending energy) of the structure. 

Before introducing the preceding assumptions explicitly, the functional in the 
form specified by both Eq. (2.1) and Eq. (2.3) can be expressed in general terms for 
a composite plate as 
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= f [n'^c® + m^K — i(n'^A*n + n^B*m+ m^B*n + m^D*m 
J-c/2[ 2 ^ 


+ 0eVl-2V:,'l) 


(2.23) 


dy — Pu^ — P77 — PwV’ — PaO 


where an integration is to be performed over the chord length of the structural 
box. The composite plate compliance matrices A*, B*, and D* are developed in 
Appendix C. Also found in the same appendix are descriptions of the resultant stress 
and moment vectors n and m, respectively. Note that the transverse shear, Vx, 
has been included in Eq. (2.23) along with a transverse shear compliance term, 0s- 
The variable 7 represents rotation of the normal to the neutral axis measured with 
respect to its initially undeformed position. In the absence of shearing deformation, 
it is assumed that normals remain normal so that 7 is equal in magnitude to the 
slope of the neutral axis. 

By introducing the assumption that the primary stresses are those that occur in 
the spanwise direction due to bending and axial stretching, the stress and moment 
resultants can be approximated as 

n = {Nx 0 0}^, m = {Mx 0 Mxyf- (2.24) 

Furthermore, considering the deformation assumptions discussed earlier in this 
section, the midplane strains and bending curvatures can be approximated by 

e® = {«o 0 0}^, = {( 7 ' + ya'') 0 2a'}^ , (2.25) 

where the prime denotes differentiation with respect to x. All of the variables in Eq. 
(2.25) are assumed to be functions of x only. Next, substituting Eqs. (2.24-2.25) 
into Eq. (2.23) and performing the indicated integration yields an expression from 
which the state vector and other terms appearing in Eq. (2.3) can be readily defined. 
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The resulting state vector is 


y = Mx Yx ^Mxy w 7 w a}^, (2.26) 

with 

— cNx, Mx — <^Mx, Y^ = cVx, Mxy = cMxy . (2.27) 

Having defined the state vector in this way, the differential equations can be deter- 
mined by the approach laid down in Section 2.1. 

It is convenient at this point to introduce a nondimension alized version of the 
resulting anisotropic beam equations. The nondimension al differential equations, in 
the form of Eq. (2.8), can be written as 

■ 0 0 0 0 0 0 0 Olf ]Vx 1 

0 0 1 00000 Mx 

0 0 0 0 0 0 0 0 Fx 

0 0 0 00000 2Mxy Pa 

An Bn 0 0 0 0 0 « 0 

Bn Bn 0 Bt3 0 0 0 0 7 0 

OOB* 00 -1 00 w 0 

Bj3 Bj3 0 B33 0 0 0 0 . , , 0 

where the nondimensional coordinate and displacement variables are 

X u w , , 

a: = -, «=£> 

It is also convenient to define 

[EI)r = crDwh , and {GJ)r = AcrD^^r, ( 2 . 30 ) 
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which are the reference values of bending and torsional stiffness in terms of the 
appropriate composite stiffnesses derived in Appendix C. (The subscript r designates 
a reference value.) In terms of the reference bending stiffness, other nondimensional 
parameters appearing in Eq. (2.28) are the force and moment resultants and external 
loads 


ct^Na^ 

M.r = 

ctMx 



{EI)jt ’ 


{EI)r ’ 

* {EI)r' 

{EI)r 

(?Pu 


^Pi 



{EI)r ’ 

Pi = 

{EI)r ' 

{EI)r ’ 

V — 

{EI)r '■ 


( 2 . 31 ) 


and the nondimensional composite compliances 


^ 11 -^^- ^11 


{EI)rD*u __{EI)RDh 

Dn = ^ , 2?13 2 ^ 


, _(EI)rBU 
2^ 

^ _ (EI)RDt^ 
^3 


( 2 . 32 ) 


and finally, the nondimensional shear compliance 


C* {EI)R^e 


( 2 . 33 ) 


It should be noted that the state vector equations presented in Eq. (2.28) 
are easily extended to both static and dynamic aeroelastic analyses by adopting 
the approach of Section 2.2 in which the load terms (i.e., inertia, damping, and 
aerodynamic) are expressed in terms of the displacement states of the structural 
state vector. The nondimensionalization of the load terms remains the same as 
that given in Eq. (2.31). Further use will be made of Eq. (2.28) when examining 


25 - 



aeroelastic behavior of composite lifting surfaces in Chapters 6 and 7. As a final 
observation about the anisotropic equations presented here, the equivalent equations 
representing isotropic structures can be obtained by simply replacing the composite 
compliance terms by isotropic compliances. 


- 26 - 


Chapter 3 


Integrating and Differentiating Matrices 


An integrating matrix approach will now be developed to solve the state 
vector equations derived in Chapter 2. The solution of such equations can often 
be a difficult task since these equations, with their boundary conditions, take the 
form of two-point boundary value problems. Additional complexity is added for 
those equations that have nonconstant coefficients. By necessity, one is forced to 
consider numerical solutions since analytical approaches can be exceedingly difficult, 
if not impossible, for all but the simplest of cases. The primary objective of this 
chapter is to discuss the development of integrating matrices, which provide the 
basis for a simple and efficient concept for numerically solving two-point boundary 
value problems. Applications of the integrating matrix to the solution of differential 
equations will be considered in the next chapter. It is hoped that the discussions here 
will lend some perspective to the integrating matrix concept as a general numerical 
tool. 

In Section 3.1, a general derivation is presented for integrating matrices that 
are suited to the integration of continuous functions. A new type of maximum 
precision integrating matrix that is developed from orthogonal polynomials will also 
be introduced. Some methods are discussed in Section 3.2 for applying continuous 
integrating matrices to the piecewise integration of discontinuous functions. Finally, 
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Section 3.3 discusses some of the aspects of a related topic, differentiating matrices. 

Although numerical solution by an integrating matrix procedure is not an en- 
tirely new idea, this approach has seen relatively little attention and development 
compared to other well known numerical tools such as finite element and finite 
difference. In fact, only a handful of investigators have contributed to the integrat- 
ing matrix method. One of the first applications of the integrating matrix method to 
problems in structural mechanics was presented in Russia by Vakhitov [26]. In this 
country, Hunter [27] is credited with much of the initial development of the integrat- 
ing matrix procedure. As mentioned by Hunter, however, an integrating matrix was 
also used by Spector [28] to simply evaluate the integrals of an asymptotic “integral 
series” solution for nonuniform beam vibration. Other major contributions to the 
application of integrating matrices (including nonlinear problems) have been made 
by White and Malatino [29] and Kvaternik, White, and Kaza [30,31]. Most of their 
analyses were for vibration and stability of rotating beams. Vakhitov [32] has also 
employed integrating matrices for a circular plate analysis, while Levashov [33-35] 
has used integrating and differentiating matrices in the context of a generalized Ritz 
method. And recently, Lakin [36] has made useful contributions to the formulation 
of integrating matrices for arbitrarily spaced grid points. 

Despite the fact that most of the applications in this work will be confined to 
structural mechanics, the developments in this chapter are quite general and can be 
applied to problems in other areas as well. It is interesting to note that the integrat- 
ing matrix technique can often be closely related to discretization methods used in 
other areas of research. To pick a single example out of many, one could reexamine 
the spline series solutions used by Schneider and Reddy [37] to solve for optimal 
nonlinear thrust vector controls for guidance of an atmospheric interceptor. This 
problem can be solved in essentially the same manner with an integrating matrix 
formulation, where the integrating matrices are developed from the appropriate 
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spline approximations. Many parallels can also be drawn with the finite element 
method. In fact, the integrating matrix approach can be considered as a special 
type of collocation finite element method. The presentation given here, however, 
will differ from the formalism of the usual finite element approach. 

As mentioned by Hunter [27], the integrating matrix also provides a useful 
tool for initial value calculations. Furthermore, the integrating matrix method 
is applicable in either linear or nonlinear situations. In the applications to be 
considered here, the focus will be on linear boundary value problems. 


3.1 Integrating Matrices for Continuous Integrands 

This section presents the basic development of the integrating matrix. As a 
preliminary requirement, it is assumed that the functions to be integrated are a 
continuous function of the spatial variable. More specifically, the integrating matrix 
development will be based upon integration of continuous polynomials that approx- 
imate the functional behavior of the structural state variables. The requirement 
of continuity, however, does not prove to be a restriction on the solution of more 
general problems. As will be shown in the next section, integrating matrices can 
be developed for piecewise continuous functions by extending the results presented 
here for continuous functions. With piecewise continuous functions, solutions are 
obtainable for almost any practical problem. 

First, a review will be given of the fundamental theory of the integrating 
matrix. This review will be independent of a specific polynomial approximation. 
The viewpoint presented here will parallel the presentations given by both Hunter 
[27] and Lakin [36], in which it is assumed that the function to be integrated 
can be represented by a polynomial of given degree. Appendix A specializes this 
fundamental presentation to the case of integrating matrices based upon Jacobi 
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polynomials. 

It should be noted that it is possible to derive many different types of integrating 
matrices, with each type dependent upon the form of approximation employed. One 
can easily formulate special purpose integrating matrices intended for a specific 
application. In this respect, the integrating matrix approach is very similar to finite 
element procedures that use element types suited to a given problem. 

The primary objective of the integrating matrix approach is to develop a 
numerical procedure for performing indefinite integrations. In contrast with initial 
value integration schemes, which are commonly used for solving differential equa- 
tions, the integrating matrix is formulated instead from a numerical quadrature. 
A quadrature is simply a numerical integration rule for integrating between fixed 
limits. An integrating matrix developed from such a quadrature rule is especially 
suited to solving boundary value problems since the region of integration is fixed 
by the boundaries. As pointed out by Hunter [27], however, the integrating matrix 
is just as easily applied to initial value problems. 

To begin, let f{x) be a continuous function on an interval [a, 6). In addition, 
suppose that a discrete set of AT -f- 1 grid points, aro, xi, ... , xn, has been chosen on 
this interval such that 


a = xo < XI < ... < XN = h, (3.1) 

and let the function values at these points be given by 

fi = f{xi). (3.2) 

In general, the points x*- can have either equal spacing or unequal spacing; ulti- 
mately, this will be determined by the nature of the approximations used for /(x). 
Furthermore, the number of points chosen is somewhat arbitrary, but for an nth 
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degree polynomial approximation to /(*) there must be at least » + 1 grid points, 
where n<N. The N + 1 grid points on [a, 6] determine N subintervals [a:,-, (t == 

0, . . . , N — 1). In general, a consecutive set of n + 1 grid points will be designated by 
the sequence Xm, i xm+n, where the subscript m denotes the starting grid 

point for the sequence. 

Assuming that f{x) can be approximated by an nth degree polynomial, then 
an approximation to f{x) on any subinterval [xi,xi^i\ can be obtained in terms of 
the values /»• = f{xi) given for a consecutive set of n + 1 grid points containing that 
subinterval. An appropriate approximation to f{x) can be obtained by any of several 
different approaches, but the most useful methods include interpolation, spline 
fitting, and least-squares fitting. In the work presented here, only the interpolation 
method will be discussed in detail. If any type of approximate data is involved, 
however, a least-squares approximation would be preferable. Lakin [36] presents 
a nice discussion of the least-squares approach as applied to the determination of 
integrating matrices. 

By integrating the approximation to f{x) over any subinterval and arranging 
the result as a linear combination of the /^’s, one obtains a convenient numerical 
description of the integration. For a typical subinterval on [a, 6] this would appear 
as 

f f{x)dx i=! + H^tm+l/m+1 H t ^tm+n/m+n] (3.3) 

Jxi 

where the If’s are weighting terms that arise from integrating the approximation to 
f{x). An approximation to f{x)dx is now easily obtained by noting that f{x)dx 
can be written as a sum of the integrations for each of the N subintervals. That is. 



+ 



f{x)dx. 


(3.4) 
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The integrating matrix representation for the approximate integration of f{x) 
contains the same information as expressed in Eqs. (3. 3-3.4), but puts it in a compact 
matrix notation suitable for matrix manipulation. First, define the column vectors 
{7} and {/} by 



and 

{/}=(/o, /l,..., fNf. (3.6) 

With these definitions, the integral in Eq. (3.3) can be expressed for all subintervals 
in the matrix notation 


(3.7) 

where the subscript n indicates the degree of the polynomial used to approximate 
f{x). “Wn is an {N + 1) X (N + 1) weighting matrix. Since the first element of {7} is 
zero, the first row of “W contains only zeros. 

A summation of the subinterval integrations can now be formally obtained 
by premultiplying both sides of Eq. (3.7) by an {N + 1) X {N + 1) lower-triangular 
summing matrix, 


$ = 


1 0 ... 01 
1 1 ... 0 

Ll 1 ... 1. 


(3.8) 


for which S,y = 1 when i>j but = 0 if t < y. As a result of this summing 
operation, if {7} is defined to be the the N +1 dimensional column vector, 

f{x) dx 

fXQ •'Xq 


( rxi rx2 rxf^f nT 

0, / f{x)dx, I f{x)dx,..., / f{x)dx\ , 

Jxn Jxn •'xn J 


(3.9) 
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then the integrating matrix relationship can be written as 

{7} = L{/} (3.10) 

where the integrating matrix L is defined by 

L=SW„. (3.11) 

Incidentally, the summation indicated by S is most easily carried out in practice by 
simple summing, rather than by matrix multiplication. 

As can be seen from Eq. (3.11), the derivation of the integrating matrix L relies 
primarily on the determination of "Wn in Eq. (3.7) because the summing matrix 
$ is known a priori. It is important to note that the integrating matrix depends 
on the polynomial approximation employed and on the number and spacing of the 
grid points, but it does not depend on the function values f{ at the grid points. 
Therefore, the integrating matrix has a separation of grid dependence and function 
dependence. Furthermore, the integrating matrix can now be viewed as a linear 
matrix operator that performs integrations via a simple matrix multiplication. 

Although the foregoing discussion has presented the general procedure to be 
followed in deriving an integrating matrix, it is now worthwhile to focus a bit more 
closely on integrating matrices that can be derived from orthogonal polynomials. 
Traditionally, orthogonal polynomials have been used as a foundation for well- 
conditioned numerical procedures. The motivation for deriving integrating matrices 
based upon orthogonal polynomials actually stems from the fact that orthogonal 
polynomials form the basis for high accuracy quadrature rules of the Gaussian 
type (cf. Section 5.4 of the text by Conte and de Boor [38]). For these Gaussian 
integration rules, the function to be integrated can be written as a product of a 
sufficiently smooth function g{x) and a nonnegative integrable weighting function 
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III 


t?(a;). That is, the integral of f{x) over (a, 6) is put into the form 

f f{x)dx= f g{x)'d{x)dx (3.12) 

J d J a 

where 

= (3.13) 

As shown in Section 5.4 of Conte and de Boor [38] and Section 7.1 of Krylov 
[39], a quadrature of the highest possible precision for a given number of grid 
points is obtained when the polynomial approximation to y(ar) is orthogonal to the 
weight function d{x) over the interval (o, 6). In addition, the nodes (grid points) are 
specified to be the zeroes of the appropriate orthogonal polynomial. These nodes 
are nonevenly spaced and are all inside of the end points o and b. 

But to be useful for the development of integrating matrices, the concept of 
an “optimal” quadrature must be extended one step further to allow nodes to be 
located exactly at the end points of an interval. The requirement for end point nodes 
becomes obvious when considering boundary value problems; these problems require 
that boundary conditions be satisfied precisely at the end points of the interval. The 
basic theory for the development of optimal quadratures having preassigned nodal 
locations is discussed in detail in Chapter 9 of Krylov [39], therefore, it will not be 
discussed in depth here. One point worth noting, however, is that fixed nodes at 
the end points give rise to a special weighting function, i?(z). This natural weighting 
function is a result of end point terms that appear in the interpolation of f{x). 

Because of the form of the weighting function that arises when end points of an 
interval are included, certain members of the Jacobi polynomial family turn out to 
be the appropriate orthogonal polynomials to use in deriving optimal quadratures. 
For this reason, the resulting integrating matrices will be referred to simply as 
“Jacobi” integrating matrices. Appendix A contains a detailed discussion of the 
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calculation procedure for Jacobi integrating matrices and Appendix B tabulates the 
corresponding weighting matrices up to n + 1 = 10. To give an idea of the accuracy 
of the quadrature on which these integrating matrices are based, all polynomials of 
degree < 2n — 1 will be integrated exactly. 

For convenient use in later calculations, weighting matrices tabulated originally 
by Hunter [27] are also repeated in Appendix B. The corresponding integrating 
matrices are referred to as “Newton” integrating matrices since they are developed 
from Newton forward difference interpolating formulas. Some additional Lagrange 
and least-squares integrating matrices not listed here are tabulated by Lakin [36], 

The Lagrange integrating matrices discussed by Lakin are somewhat related 
to the Jacobi matrices described above since both originate from Lagrange inter- 
polations, which are valid for unequal grid point intervals. As noted by Lakin, 
however, the Lagrange integrating matrices are somewhat cumbersome to numeri- 
cally compute for grid spacings chosen on an ad hoc basis. In contrast, however, 
the computation of Jacobi matrices is a much simpler numerical task. It turns out 
that Jacobi integrating matrices can be calculated by a procedure that is in many 
respects similar to the procedure presented by Lakin for least-squares fitting based 
on orthogonal polynomials. In fact, if the least-squares fit procedure is applied to 
the Jacobi grid points for an approximating polynomial of the maximum degree (i.e., 
the same degree used for an interpolation), then the least-squares procedure yields 
the Jacobi integrating matrix. A nice feature of the Jacobi integrating matrices is 
that optimal grid point locations are determined automatically by the underlying 
quadrature rule. 

There are some other features of Jacobi integrating matrices that differ from 
Newton and Lagrange integrating matrices. The first of these differences arises 
because of the unequal grid point spacing. In the calculation of Newton matrices 
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for instance, the interpolations over a subinterval [x^, are performed by roving 
polynomials. That is, a polynomial of degree n, which makes use of n+1 consecutive 
grid points, can be shifted along the interval of iV + 1 grid points, one point at 
a time. This is possible, of course, because the Newton matrices are based on 
equal grid point spacings. For the interpolations required by the Jacobi integrating 
matrices, however, the unequal subinterval lengths mandate the use of a stationary 
polynomial. Because of the use of stationary polynomials, the Jacobi integrating 
matrices are in some respects analogous to high order, polynomial-based finite 
elements. That is, a Jacobi “element” corresponds to n + 1 consecutive grid points 
and has n — 1 unequally spaced internal nodes. The complete interval [a, 6] can be 
constructed by placing so-called Jacobi “elements” end to end. Experimentation 
with the Jacobi matrices reveals that the highest numerical efficiency is obtained 
by using a small number of “elements” of high order. This is in accord with results 
for finite element and other numerical approximation techniques. 

A second aspect of a Jacobi integrating matrix (or any other type of integrating 
matrix with unequal intervals) is that interpolation may be required if one desires 
solution results at points other than the grid points. This is a fairly simple process, 
however, since interpolation shape functions are easily developed for the Jacobi 
polynomials. These shape functions are presented in Appendix A. Again, these 
shape functions are analogous to shape functions commonly used in finite element 
analysis. 

Regardless of the type of integrating matrix, some very useful information can 
be obtained by investigating the quadrature rule on which the integrating matrix 
is based. The quadrature rule consists of the weighting terms that are applied to 
the integrand at each of the grid points. To be specific, the quadrature rule for the 
integral 
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simply consists of the W’s and the corresponding a:,’s. The W’s, which are the 
primary values of interest, are contained in the last row of the integrating matrix. 
These quadrature weights, of course, can also be obtained by summing each of the 
terms in a column of the weighting matrix W. 

An examination of the underlying quadrature can often lead to a better un- 
derstanding of the nature of a particular integrating matrix. For example, for the 
Jacobi weighting matrices given in Appendix B, the grid points are arranged sym- 
metrically on the normalized interval (—1, 1] and the quadrature weights (obtained 
by summing the columns) are all positive numbers that steadily increase in value 
as one approaches the midpoint of the interval. The values of the weight terms are 
also symmetric about the midpoint of the interval. But the important revelation is 
that the quadrature rule associated with Jacobi integrating matrices turns out to 
be the same as a well known numerical integration method, Lobatto quadrature. 
In fact, the numerical calculations for the Jacobi weighting matrices in Appendix B 
were verified by comparing the quadrature weights with those tabulated for Lobatto 
integration on page 920 of Abramowitz and Stegun [40]. 

An aspect of Newton integrating matrices that became apparent during this 
investigation was a possible asymmetry of the underlying quadrature. This phenom- 
enon, which was also noted by Lakin [36], is referred to as “biasing.” Biasing only 
arises when using interpolating polynomials for which the number of interpolation 
points, n -t- 1, is an odd number. The reason that biasing occurs for odd n -f 1 is 
that away from the end points the interpolating polynomials cannot be centered 
on the subintervals [a:t',*t+i] that are being interpolated. When biasing is present, 
integrating matrices corresponding to Newton forward difference formulas differ 
from those derived from backward difference formulas. When n -i- 1 is even, there is 


no biasing and both forward difference and backward difference formulations yield 
identical integrating matrices. When using the Newton integrating matrices, it 
is generally convenient to employ those matrices that are not biased (i.e., based 
upon TP„ with n + 1 even). It is interesting to note, however, that there is an easy 
way to symmetrize a biased Newton integrating matrix. Essentially, one develops 
a weighting matrix from one-half the sum of forward difference and backward 
difference weighting matrices. This is conveniently illustrated by taking a particular 
example that considers a quadratic Newton weighting matrix written for five grid 
points. By expressing the weighting matrix as a sum of forward and backward 
difference matrices, we have 


'^2 = ^(W2,/ + W2,&) 


( 3 . 15 ) 


where the subscripts / and b refer, respectively, to forward and backward differ- 
encing. Equation (3.15) expands into (see Appendix B) 



'0 

0 

0 

0 

O' 


■ 0 

0 

0 

0 

or 


5 

8 

-1 

0 

0 


5 

8 

-1 

0 

0 


0 

0 

5 

0 

8 

5 

-1 

8 

0 

-1 

- 1 - 

-1 

0 

8 

-1 

5 

8 

0 

5 

0 

0 


0 

0 

-1 

8 

5 


0 

0 

-1 

8 

5 


.5 

13 

11 

15 

4J 


L 4 

15 

11 

13 

5 J , 


( 3 . 16 ) 


■ 0 

0 

0 

0 

O ' 

10 

16 

-2 

0 

0 

-1 

13 

13 

-1 

0 

0 

-1 

13 

13 

-1 

0 

0 

-2 

16 

10 

L 9 

28 

22 

28 

oJ 



where the values below the horizontal bar in each matrix represent the quadrature 
weights. As can be seen, the final weighting matrix in Eq. (3.16) yields a symmetric 
quadrature. An interesting comparison can be made with the equivalent nonbiased 
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cubic Newton weighting matrix. It appears as 


W3 
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( 3 . 17 ) 


The similarity of the interior matrix elements between Eq. (3.17) and the last expres- 
sion in Eq. (3.16) indicates that the symmetrized, nonbiased quadratic matrix will 
approach the accuracy of the cubic matrix. Obviously, error cancellation inherent 
in the summing of forward and backward difference formulations is responsible for 
the accuracy increase. 

Also, rough comparisons between integrating matrices can be made based 
upon the quadrature rules. For instance, the Newton integrating matrices, in 
contrast with Jacobi integrating matrices, have quadrature weights that can oscillate 
considerably as one proceeds along the interval of integration. This oscillation tends 
to increase somewhat for the higher order Newton integrating matrices. This basic 
difference between the underlying quadratures for Newton and Jacobi integrating 
matrices can manifest itself in both accuracy and convergence properties. This will 
be discussed again in later chapters when comparisons are made between Newton 
and Jacobi solutions. It should be noted, however, that Newton and Jacobi matrices 
are identical for the quadratic approximation. For higher order approximations, 
their properties and relative accuracies differ. Although both types of integrating 
matrices provide very good solution accuracies, it has been found that Jacobi 
matrices are capable of offering faster, more predictable convergence and higher 
accuracy with fewer grid points. 
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3.2 Integrating Matrices for Discontinuous Integrands 

Aided by the developments of the previous section, it is now straight forward 
to determine appropriate weighting matrices and integrating matrices for systems 
with discontinuous parameters. Basically, there are two types of discontinuities that 
occur often in static structural analysis and structural vibration. The first of these 
discontinuities stems from stepwise changes in coefficient terms of the differential 
equations. Such stepwise discontinuities could easily be the result of changes in 
stiffness or mass parameters. Appropriately, these discontinuities are best described 
by step functions. Methods for handling this type of discontinuity will be considered 
in this section. The methods discussed here for treating stepwise discontinuous 
systems are similar in some respects to the methods described by Vakhitov [26]. 

A second type of discontinuity results from point loading. Point loads can arise 
from either applied external loads or inertia loads associated with point masses. 
Point loads are most easily handled by introducing impulse functions. Since impulse 
functions are to be treated by an approach different from that for step functions, the 
methods for including concentrated loads in an analysis are taken up in Appendix 
E rather than this section. 

The treatment of a problem with stepwise discontinuities is straightforward. 
The solution interval is simply divided up into analytic regions by breaking the 
integration at points of discontinuity. For each of these separate regions a different 
(or the same) weighting matrix can be used. An integrating matrix for the complete 
solution interval is then obtained in the normal way by summing the weighting 
matrix terms for each of the piecewise regions. This summation procedure is 
identical with the one carried out by the summing matrix S in Eq. (3.11). The 
only difference in forming the integrating matrix for stepwise discontinuous systems 
is that the weighting matrix “W for the complete solution interval has a block matrix 
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Fig. 2. Typical integration regions for a discontinuous function 


structure, where each of the blocks corresponds to a normal weighting matrix 
written for one of the analytic regions. 

Perhaps a simple example best explains the procedure to be applied to dis- 
continuous functions. Consider integrating a function that has a single point of 
discontinuity, as shown in Fig. 2. The integration of the discontinuous function can 
clearly be broken into the two distinct integrations indicated by the regions I and 
II. For this particular example, which assumes weighting matrices with end point 
nodes, the equivalent of Eq. (3.7) can be written as 




TV/ 

. 0 


0 ]{/,| 


= "Wi+Iif. 


( 3 . 18 ) 


Here, as in Eq. (3.11), the resulting integrating matrix can be calculated from 


L = SWi+n. ( 3 . 19 ) 

In Eq. (3.18), it should be noted that the “merged” weighting matrix has no overlap 
of rows or columns between the individual matrices for each region. 


- 41 - 


Note in Fig. 2 that two distinct function values have been defined at the discon- 
tinuous point, each corresponding to the appropriate function value an infinitesimal 
distance away from the discontinuity. These identically located points occur for the 
case in which the weighting matrices on either side of the discontinuity make use 
of grid points located at the end points of an interval. Strictly speaking, it is not 
mandatory that one have grid points at the discontinuity, unless, of course, bound- 
ary conditions or loads must be applied at such a point. For instance, an interior 
region of a multiple region integration might make use of weighting matrices based 
on Gauss-Legendre quadrature, which does not require grid points at the ends of 
an interval (see Appendix B). Regardless of whether or not grid points are located 
at discontinuities in the integrand, the continuity of the integration assures that 
solutions will be continuous. 

Another feature that can be accommodated in much the same fashion as a 
stepwise discontinuity in integration is either a change in step size or a change in 
integrating matrix type within a solution interval. Grid points are fixed, of course, 
for weighting matrices based on orthogonal polynomials, but can be variable for 
constant step size matrices like the Newton integrating matrices. It is quite easy 
to switch to different types of integrating matrices to satisfy particular solution 
requirements for portions of a solution interval. The method for handling a discon- 
tinuous integration arising from integrating matrix changes or step size changes is 
very much similar to the method employed for discontinuous system parameters. 

When compared to the discontinuous problem in Eq. (3.18), changes in either 
step size or integrating matrix type lead to a slightly different way of adding the 
local weighting matrices into the global weighting matrix. If weighting matrices 
with different step sizes (or even with the same step sizes) are merged to integrate a 
larger region, there will be some “overlap” between weighting values for the merged 
matrices. This is true for weighting matrices that use end point nodes as well as 
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Fig. 3. Typical regions for piecewise integration of a continuous 
function 

those that do not. The “overlap” shows up in the global weighting matrix as an 
overlapping of some of the rows and columns of adjoining matrix blocks. Again, 
a specific example easily demonstrates the procedure. Consider the five grid point, 
two region integration shown in Fig. 3. For the sake of illustration, second degree 
Newton weighting matrices having equal step sizes will be employed in both regions 
I and II. Furthermore, the the step size is conveniently chosen such that the 
constant step size factor h/l2 has the value of unity. (See Appendix B.) With these 
specifications, the equivalent of Eq. (3.7) becomes 
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(3.20) 


For this specific example, it can be seen that the region I and region II submatrices 
overlap by one row and one column. 
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In the integrating matrix method, the process of merging weighting matrices 
is analogous to the merging of stiffness matrices in the finite element method. In 
fact, by considering the different situations that can arise in connection with discon- 
tinuous integrations, it is possible to write down a simple set of rules describing how 
weighting matrices should be merged. To be able to give a consistent set of rules for 
weighting matrices with and without end point nodes, one basic requirement must 
be met. This requirement simply stipulates that weighting matrices without end 
point nodes must conform to weighting matrices with end point nodes. An example 
of a conforming Gauss-Legendre weighting matrix is given in Appendix B. It should 
be noted too that the merging process is carried out by starting at one boundary 
of the solution interval and proceeding with consecutive piecewise regions in the 
direction of integration. 

The rules for merging weighting matrices at points of continuous or discon- 
tinuous system coefficients can be considered separately. First, for discontinuous 
parameter systems we have the following requirements for merging at a point of 
discontinuity: 

1. If both adjacent matrices have end point nodes, there will be 
no row or column overlap (e.g., Eq. (3.18) ). 

2. If one of the adjacent matrices has end point nodes and the 
other does not, there will be one row and one column of 
overlap. 

3. If neither of the adjacent matrices have end point nodes, there 
will be one row of overlap and two columns of overlap. 

For continuous parameter systems, the following rules for merging apply at any 
point within the solution interval: 

1. If either or both of the adjacent matrices have end point nodes, 
there will be one row and one column of overlap (e.g., Eq. 

(3.20) ). 

2. If neither of the adjacent matrices have end points, there will 
be two rows and two columns of overlap. 
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Any of the above merging rules can be easily checked by taking a sample 
problem and examining the subinterval integrations, as demonstrated, for example, 
in Eq. (3.20). Since the merged weighting matrices only specify the subinterval 
integrations, the logical conclusion of the merging process is in the determination 
of the corresponding integrating matrix from Eq. (3.11). From the standpoint of 
numerical calculation, the merging rules can be applied during the summing process 
that forms the integrating matrix from the weighting matrix. 


3.3 Differentiating Matrices 

Similar to the idea of an integrating matrix, there is also a differentiating 
matrix. The differentiating matrix expresses the derivative of a function in terms 
of the function values at a discrete set of points. Usually, a differentiating matrix 
is determined by differentiating an interpolating polynomial, just as the integrating 
matrix is determined by integrating that same polynomial. It is also feasible to use 
simple differencing to calculate a differentiating matrix. 

Differentiating matrices, in fact, are quite common in finite element applica- 
tions. For instance, the strain- displacement relationships in a finite element analysis 
are usually written in terms of a differentiating matrix determined by differentiating 
the displacement shape functions. The papers by Levashov [33-35] show how a 
combination of both differentiating matrix and integrating matrix concepts can be 
used to solve structural problems by the displacement method. For the integrating 
matrix approach as presented here, one rarely has to resort to using differentiating 
matrices. Nonetheless, the differentiating matrix proves to be a useful tool in certain 
circumstances and is worthy of mention. 

With the state vector formulation of the aeroelastic equations given in Chapter 
2, it is normally possible to write the structural portion of the equations without 
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need for a differentiating matrix. To express certain descriptions of the aerodynamic 
loads, however, higher derivatives of some of the structural state variables may be 
required. In such a situation, the differentiating matrix provides a convenient way 
of writing these higher derivatives in terms of existing state variables. A similar 
situation can exist when writing the “preload” terms for a structure with initial 
curvature. In this latter case though, it appears that it’s usually possible to rewrite 
the equations in such a way that the use of a differentiating matrix is avoided. 

A variety of differentiating matrices can be calculated depending on the type of 
interpolating polynomial assumed. Methods for calculating differentiating matrices 
are quite straight forward since one simply determines derivative expressions for 
the approximating polynomials. Detailed calculations will not be presented here for 
differentiating matrices; however, some precalculated formulas do exist for deriva- 
tives of certain types of interpolations. For example, expressions for the derivatives 
of Newton forward difference interpolation formulas, which form the basis for the 
Newton integrating matrices, can be found on page 883 of Abramowitz and Stegun 
[40]. Also tabulated there are similar expressions for the derivatives of Lagrange’s 
interpolation formula, which is useful for approximating functions with nonuniform 
grid spacings. And finally, differentiating matrices corresponding to Jacobi integrat- 
ing matrices can be calculated by differentiating the shape functions presented in 
Appendix A. The finite series expansion for the Jacobi polynomials presented in 
that appendix should prove useful in obtaining the proper derivative expressions. 
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Chapter 4 


Integrating Matrix Solution of 
State Vector Equations 


Integrating matrices provide a convenient and efficient method for 
solving two-point boundary value problems. The prime concern is the solution of 
the aeroelastic and structural dynamic equations presented in Chapter 2. The major 
focus of this chapter is on the formulation of discrete solutions for these equations 
when they are represented in state vector form. Other important issues that are 
discussed include methods for manipulating and reducing the discrete state vector 
equations to a suitable form for numerical solution. 

A direct approach that does not require the calculation of influence matrices will 
be presented for formulating linear systems of equations and eigenvalue problems. 
It also turns out that the integrating matrix solutions to be discussed are closely 
related to alternate solutions in terms of transition matrices. In fact, the integrating 
matrix provides one of the easiest means for numerically calculating transition 
matrices. Furthermore, with the aid of a transition matrix, one can in turn calculate 
symmetric stiffness influence matrices for the types of problems considered in this 
chapter. Although the transition matrix approach will not be discussed here, it is 
presented briefly in Appendix G. Many of the ideas contained in this chapter can 
be used in connection with transition matrix calculations. 



As will be evident to those familiar with approximate solution techniques, the 
integrating matrix approach belongs to a class of collocation solution methods. 
Since collocation solutions belong to the more general family of methods know as 
weighted residual methods, it is possible, in fact, to formalize the integrating matrix 
approach in much the same way that finite element methods are formalized. For the 
sake of simplicity, however, the following presentations will employ a more intuitive 
approach. 

Collocation solutions are conveniently classified by Collatz [41] (also, Meirovitch 
[42]) according to the nature of the approximating functions. Within this classifi- 
cation scheme, the integrating matrix approach presented here falls into the category 
of a mixed method. That is, the approximating functions initially satisfy neither the 
differential equations nor the boundary conditions of the problem. (Other categories 
stipulate that the functions must initially satisfy either the boundary conditions or 
the differential equations.) Solutions with the mixed method are obtained by insist- 
ing that at selected collocation grid points the boundary conditions, or diflferential 
equations, as appropriate, are satisfied exactly. From the discussions of chapter 3, 
it should be clear that the arrangement of collocation grid points, as well as the 
approximating functions, are specified by the choice of integrating matrix. 

A general advantage that collocation methods have over some other methods 
lies in the ease with which coefficient matrices can be generated for the approximate 
equations. As will be shown, the integrating matrix approach provides a very simple 
way of obtaining the matrix equations that describe a particular set of differential 
equations. In addition, if these differential equations are in state vector form, it is 
often possible, at least for linear problems, to take advantage of the structure of 
the state vector equations to further simplify the corresponding matrix equations. 

There is, however, a disadvantage to collocation type approaches. Namely, 
the matrices involved are always nonsymmetric, thereby requiring less efficient 
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solution routines. But in many instances, this disadvantage is entirely offset by 
the simplicity of the approach. Furthermore, since the integrating matrix method 
easily accommodates very high order approximating functions, one often finds a 
considerable reduction in problem size for a required solution accuracy. 


4.1 Discretized State Vector Equations 

To obtain a discretized form of the state vector equations, it is convenient to 
consider a particular example problem. A suitable choice is Eq. (2.20), which can 
be used to describe aeroelastic response for a one dimensional structure. Although 
this equation is intended primarily for static aeroelastic response, with appropriate 
generalization it encompasses other types of equations, including dynamic equations 
like Eq. (2.21) and Eq. (2.22). To make the results generally applicable, it is assumed 
that Eq. (2.20) can be expressed in a nondimensional form so that it appears as 

y' = Zy-XAy-Sr (4.1) 

where X represents a nondimensional parameter. For homogeneous problems, which 
are obtained by dropping the nonhomogeneous term ar, the parameter X is to be 
taken as an eigenvalue. In the sample problems to be presented later, more precise 
definitions will be given for the terms appearing in Eq. (4.1), but for now it is 
convenient to use them simply as generic terms. In what follows, the state vector y 
in Eq. (4.1) will be referred to as a local state vector. 

Since Eq. (4.1) must be valid for any value of the spatial coordinate, it can be 
written at a discrete set of iV + 1 points along the desired solution interval. In this 
discrete form it can be expressed as 

y' = Zy - XAy - ir- (4.2) 
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The discrete version of the state vector, y, will be referred to as a global state 
vector. Note that the dimension of y will be NS(AT + 1), where NS refers to the 
number of states in the local state vector. The ordering of components in the global 
state vector is actually arbitrary, but it turns out to be very useful in what follows 
to specify a particular ordering. It will be assumed that the global state vector 
is ordered such that its structure is the same as that of y. More precisely, it is 
partitioned into generalized force and generalized displacement subsets as was done 
for the local vector y in Eq. (2.4), and furthermore, the discrete set of values for a 
particular state variable are grouped together. This will be clearly demonstrated in 
the sample problems of Chapter 5. 

With the aforementioned ordering in mind, a global version of the integrating 
matrix, L, can now be applied as a matrix operator to both sides of Eq. (4.2). The 
result is 


y = LZy — XLAy — Lir + k (4.3) 

where k is a constant vector of integration to be determined from boundary con- 
ditions. Note that the effect of applying the integrating matrix to the isolated 
derivative term on the left-hand side of Eq. (4.2) was simply to integrate that term 
and thus remove the derivative. Integrating matrices on the right-hand side of Eq. 
(4.3) perform their integration function through a matrix multiplication. 

The global integrating matrix in Eq. (4.3) is a block diagonal matrix that 
appears as 


L = 


(4.4) 


with each block being a standard (N -i- 1) X (N 1) integrating matrix as derived 
in the previous chapter. The appropriate number of matrix blocks on the diagonal 


- 50 - 



is determined by the number of state variables, NS, in the local state vector. In 
essence, the integrating matrix allows one to proceed in much the same manner as 
solving a simple integrable equation by analytical methods. But with the numerical 
integration properties of the integrating matrix, it is now possible to deal with 
otherwise difficult integrations in a very simple manner. 

At this point, Eq. (4.3) provides a matrix equation from which one can obtain 
a solution of the state vector differential equations. The only remaining step is 
to determine k from specified boundary conditions. And once k is known, there 
remains only a linear systems solution if considering the nonhomogeneous problem, 
or a matrix eigenvalue solution if considering the related homogeneous problem. 

For two-point boundary value problems, it is useful to introduce the boundary 
condition matrices Bq and Bn. These boundary condition matrices, which relate 
to homogeneous boundary conditions, will aid in solving for k. Bq and B„ can be 
written as 

Bo = Ibg’ (4.5) 

and 

Bn = ibl (4.6) 

where 1 is a column vector containing all unit terms. The vectors bo and b„ can be 
expressed as 


bj = {l,0,...,0) 


(4.7) 


and 


bj[ = {0,...,0,l}. 


(4.8) 
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The dimensions of Bq and are (iV+ 1) X {N + 1), the same as a normal integrating 
matrix. Also, since the first row of the integrating matrix contains only zeroes, one 
has that BqL = 0. 

Corresponding to the global block diagonal integrating matrix L, there is a 
similar block diagonal boundary condition matrix, B. For a two-point boundary 
value problem, each matrix block on the diagonal of B is specified by applying either 
Bo or Bn (or equivalently, bj or b^) to the corresponding state variable and solving 
for its constant vector of integration. This process will be clearly demonstrated in 
the examples of Chapter 5. Since bg and b^ contain mostly zero elements, along 
with a strategically located unit term, their operational effect on a discrete state 
vector is to select the “degree of freedom” at which a boundary condition is to be 
applied. 

In addition to the homogeneous boundary condition matrix B, one can define 
B„/i to account for nonhomogeneous boundary conditions that can be written in 
terms of the state variables. The specific form of B„/, has to be determined for each 
particular problem, but one should note that for many common problems it is simply 
zero. When B„/i does need to be determined, it is defined in such a way that when 
it premultiplies the global state vector, it produces the required nonhomogeneous 
boundary terms. Similar to the situation for B, B„/, will consist mainly of zeroes, 
but will have a few strategically located nonzero terms. The nonzero terms in B„;,, 
however, are not usually unit terms as was required for B, 

With the foregoing definition of boundary condition matrices, a general expres- 
sion can be obtained for the solution of k. To obtain this expression, Eq. (4.3) is 
first multiplied through by B. Since B has been defined for homogeneous boundary 
conditions, we have that By — 0. Furthermore, the form of k is specified to be such 
that Bk = k. With these two identities, and with the aid of the nonhomogeneous 
boundary term B„;iy, one obtains from Eq. (4.3) the general result 
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k = -BL(Zy - XAy - ir) - B„/,y 


(4.9) 


where it can be noted that k consists of both homogeneous and nonhomogeneous 
components. For nonhomogeneous boundary terms that express displacement de- 
pendent loads, it is often possible to rewrite these terms is such a way that they can 
be included in A rather than B^h- The advantage in doing this will be discussed in 
the next section. 

Having now determined a general expression for k, the next step is to obtain 
a suitable matrix equation that can be solved by standard methods. This is easily 
achieved by substituting the expression for k from Eq. (4.9) into Eq. (4.3). Grouping 
similar terms and then rearranging leads to 

[H-XFA]y = f (4.10) 

where 


H = I-hB„;,-tFZ 

(4.11) 

F = [B - IJL 

(4.12) 

f = Fir , 

(4.13) 


with 1 being an identity matrix with appropriate dimensions. It is obvious that 
Eq. (4.10) represents a system of linear equations when X is specified and when the 
nonhomogeneous external load term, f, is nonzero. If f is zero, then Eq. (4.10) 
provides a matrix eigenvalue problem. Furthermore, if one chooses to give different 
interpretations to XA, then Eq. (4.10) encompasses a broad range of aeroelastic, 
vibration, and structural problems. For instance, the free vibration and unsteady 
aeroelastic problems described by Eqs. (2.21-2.22) are manipulated into similar 
forms after application of the integrating matrix. More precise definitions for several 
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types of problems will be given in the following sections and in the example solutions 
of Chapters 5, 6, and 7. 

It is possible to stop at this point and numerically solve Eq. (4.10) with standard 
methods. To do so, however, requires the solution of matrix equations with twice the 
number of necessary variables. As mentioned previously, the global state vector, y, 
is partitioned into both force and displacement variables, and y^,, similar to Eq. 
(2.4). The next section will show that by applying matrix partitioning techniques 
one can conveniently reduce the solution of Eq. (4.10) to an equivalent problem 
written solely in terms of the displacement variables y^,. 


4.2 State Vector Equation Reduction 


By partitioning the matrices in Eq. (4.10) according to the force and displace- 
ment states of the global state vector, it is possible to achieve an analytical reduction 
of the matrix equations such that only displacement variables are involved in the 
final solution. As will be shown, the reduction task is greatly simplified due to the 
structure of the state vector equations developed in Section 2.1. 

Writing out the partitioned form of the matrices appearing in Eq. (4.10) reveals 

that 


H = 


Hff 

y^DF 


^FD 

^DD. 


(4.14) 


A = 


0 

.0 


^FD 


0 J 



(4.15) 


(4.16) 



and 



First of all, it should be noted that the H matrix contains structural related 
terms, including elastic restraint boundary terms that arise from B„/i. The Wfd 
submatrix consists of these elastic restraint terms as well as “geometric stiffness” 
terms for problems involving initial structural curvature and deformation dependent 
“preloads.” As mentioned in connection with nonhomogeneous boundary terms in 
the previous section, it is often possible to shift terms that appear in Wfd so that 
instead they can be included in Afd- This can be easily done, for example, for the 
centrifugal stiffening term that appears in rotating beam problems. The motivation 
for rearranging in this way is to zero out the submatrix. If Wfd = 0, the 
reduction process is considerably simplified, as will be indicated in what follows. 
For problems without nonhomogeneous boundary terms and without geometric 
nonlinearities, Wfd is automatically zero. 

Other matrices and vectors appearing in Eqs. (4.15-4.17), such as A and f, have 
the indicated forms because of the natural structure of the state vector equations 
derived in Section 2.1. For instance, since A contains only displacement dependent 
load terms, the only nonzero partition is given by ^fd, which is the submatrix that 
expresses loads in terms of the displacement variables. F, on the other hand, is a 
block diagonal matrix simply because the matrices from which it is calculated are 
also block diagonal. 

With the partitioned form of the matrices given by Eqs. (4.14-4.17), it is 
straight forward to write out the corresponding partitioned equations from Eq. (4.10) 
and use them to eliminate the the generalized force variables, y^. This elimination 
process yields the expressions to be used in solving for the force variables once 
the displacement variables are obtained. The method to be used in reducing 
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the partitioned equations involves multiplying them through by the inverse of the 
H matrix, where this inverse is obtained from the partitioned form of H. Matrix 
inversion by partitioning is discussed, for example, in Appendix A of the text by 
Przemieniecki [43]. If the inverse of H is denoted by 






FF 




H * 1 

FD 


H* 


DD-i 


(4.18) 


then the inverted terms that will be needed in the following reductions are 


h;^/. = h 


-1 

FF 


+ \^DF^fF 


(4.19) 


and 


X-1. 


= -{»DD - HzjfHJ/. 


(4.20) 


The first set of reduced equations to be obtained are those for aeroelastic lift dis- 
tribution and structural deflection problems. Beginning with the nonhomogeneous 
linear equations given in Eq. (4.10), their partitioned form appears as 


”Hff 

Hfd' 

- X 

0 FFFAFDll|yF| 


l^DF 

^DD. 


o 

0 

1 

1 0 1 


(4.21) 


Multiplying Eq. (4.21) through by H ^ then yields 


I 

.0 


I — lyi)/ 


^FF»rp 


(4.22) 


By applying Eqs. (4.19-4.20), this can be simplified to give the generalized force 
equation 

Vf — — — XF^FAFoIXi, + FF^rf (4-23) 
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and the reduced linear system 


I — XTAf£> ~ 


Vd 


Ti 


rp 


(4.24) 


where 


T = 


~^DD^DF^FF^ ff ■ 


If Hfx» = 0, then Eq. (4.24) simplifies to 


(4.25) 


[I XTAfx>]y£) — . (4.26) 

This simplified result provides the expression to be used in several examples that 
will be presented in later chapters. Note that X must be specified before solutions 
can be obtained from Eqs. (4.24) and (4.26). In the absence of aerodynamic or 
inertial loading, X will be equal to zero and Eq. (4.26) directly provides structural 
deflection solutions. 

The next important result concerns the reduced eigenvalue equations for both 
aeroelastic divergence and free vibration problems. This result is easily reached by 
first setting the right-hand side of Eq. (4.21) to zero to obtain the homogeneous 
version of the partitioned equations. Then, by following an approach similar to 
that given above for the nonhomogeneous case, one obtains the reduced eigenvalue 
equations 


AiP£> — (1/X)l]y£, = 0 . 


(4.27) 


The corresponding generalized forces, after solving for y^,, are given by 


ip = XHi^^F^jT-AF/jy^ . 


(4.28) 
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Similar to Eq. (4.26), a useful form of the reduced eigenvalue problem in Eq. (4.27) 
is obtained by considering = 0, which gives the result 

[TAf£> — (1 A)I]y£) = 0 . (4.29) 

This form of the eigenvalue problem will be employed in later examples. 

To complete the presentation of the reduced equations, it is useful for reference 
purposes to give the flutter eigenvalue problem that will be discussed in Chapter 7. 
The flutter equations, which will be written for the Laplace domain, are obtained 
in a manner similar to Eq. (4.29), except that they derive initially from Eq. (2.21) 
rather than Eq. (2.20). By taking Hfd = 0, these equations can be written in the 
form 

I + + CpDS* — Qfz?(s*, X)] y/) = 0 (4.30) 

where is a mass matrix, Cfd is a damping matrix, and is an unsteady 
aerodynamic matrix. The variables s* and y^, are, respectively, a nondimensional 
Laplace variable and a nondimensional complex eigenvector. The terms appearing 
in Eq. (4.30) will be defined in more detail in Chapter 7. 

An important observation about the foregoing reductions is that although the 
required matrix inversions can be carried out numerically, they can often easily 
be avoided or simplified. This is possible because of the particularly simple block 
matrix structure of H. For example, if Hfd = 0, then the inversion of H is given by 

Hff 0 I"* r Hp/ 0 

.Hx>f HuJ 

Furthermore, as a result of the structure of the hybrid state vector equations, 
inversions for Wdd and H/?/? are trivially obtained. This is due in part to the 


( 4 . 31 ) 
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sparseness of these matrices. Further clarification of the inversion process will be 
offered in the examples of Chapter 5. As a final result, however, one will find it 
possible to write out the simple matrix expressions for the calculation of T and its 
product with other matrices. 

Additional simplification of the reduced equations is possible for situations 
where constraint equations exist. The equations of constraint that will be applied 
in later chapters express simple relationships between solution variables. Such 
relationships between variables can be used to eliminate certain degrees of freedom 
by writing them in terms of other degrees of freedom. The methods for applying 
such constraints to the discrete state vector equations are presented in Appendix F. 
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Chapter 5 


Solutions for Isotropic Beams and Rods 


Integrating matrix solutions are most easily demonstrated with the 
aid of sample problems. This chapter, and the next two, present examples of 
how integrating matrix solutions can be used to solve a variety of problems in 
structural dynamics and aeroelasticity. These solutions are presented both to 
validate the methods discussed in Chapter 4, and to evaluate the accuracy and 
convergence trends of the integrating matrix solutions. Comparisons are made 
between the integrating matrix solutions and alternate analytical or numerical 
calculations. Among the items considered in the following examples are various 
types of boundary conditions and loading conditions, all of which commonly arise 
in performing practical vibration studies and aeroelastic analyses. 

5.1 Axial Vibration of Cantilevered Rods 

Axial vibration of a rod provides a simple example with which to study the 
integrating matrix solution process. Since axial rod vibration is described by a 
second order Sturm-Liouville differential equation, analytical solutions are available 
for comparison with numerical solutions. The integrating matrix solution properties 
presented for these simple Sturm-Liouville systems carry through to the much more 
complicated problems to be discussed later. Aeroelastic divergence of unswept 
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isotropic wings is also described by a similar Sturm-Liouville differential equation. 
Therefore, the fundamental eigenvalue of the following eigenvalue solutions can 
alternately be interpreted in terms of divergence dynamic pressure. 

A nondimensional form of the second-order Sturm-Liouville problem describing 
axial rod vibration can be written as 

d { \ du\ ^ , . 

where 




r(s) = 


{EA)n 
EA ’ 


oj^rriRt^ 

{ea)r • 


niR 


(5.2) 


For this simple differential equation, the dimensionless state vector form can be 
written down directly, yielding 


where 



(5.3) 


F«(ir) = 


1 du 

T d^ ' 


(5.4) 


If the rod is considered to be clamped at the end x = 0, the cantilever boundary 
conditions for Eqs. (5.1) and (5.3) can also be specified as 


u(0) = 0, F„(1) = 0. (5.5) 

The differential equations appearing in either Eq. (5.1) or Eq. (5.3) can be readily 
solved with an integrating matrix, as will be demonstrated next. 
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5.1.1 Rods with Continuous Mass and Stiffness 


The first solution to be presented will consider axial vibration for a cantilevered 
rod that can have variable mass and stiffness parameters. Because of its simplicity, 
this problem serves as a good demonstration of the integrating matrix solution 
procedures. The integrating matrix will be applied first to the hybrid state vector 
equations given in Eq. (5.3). It is also instructive to demonstrate how the integrating 
matrix can be utilized to solve the differential equations given in the form of Eq. 
(5.1). Both types of solutions will be presented in detail to serve as models for other 
problems discussed in this chapter. 

Beginning with the discretized form of the state vector equations in Eq. (5.3), 
the integrating matrix is applied as an operator to perform the necessary integra- 
tions. The next step is the application of boundary conditions to determine the 
boundary condition matrices B and B„/,. From the boundary condition matrices 
follow the definitions for the other matrices appearing in Eqs. (4.10-4.13), which 
have the partitioned forms displayed in Eqs. (4.14-4.17). For example, in a homoge- 
neous vibration problem, for which ir = 0, the F matrix is defined in terms of the 
given boundary condition matrices, followed by the determination of the H matrix. 

For the axial vibration problem, an application of the global integrating matrix 
to the discretized version of Eq. (5.3) yields 


D 'l o-ifro fo 

.0 L. I ,'r 0. Iw/ .0 O.ltiJJ lk„j 


( 5 . 6 ) 


To make it convenient to solve for boundary conditions, Eq. (5.6) is easily expanded 
into the two sets of equations 


{Fu} = —XL' m{«) -t- kf 
{«} = L'r{F„}+ k„. 
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( 5 . 7 ) 

( 5 . 8 ) 


Applying the boundary condition «(0) = 0 to Eq. (5.8) gives 


bo{a}=0 = ku ( 5 . 9 ) 

where use has been made of b^’L = 0 and bj’ku = A:„. (Note, for example, that 
l(„ = iku, where 1 is as defined in Eqs. (4.5) and (4.6).) In a similar manner, the 
boundary condition F„(l) = 0 can be applied to Eq. (5.7), yielding the scalar solution 

bn{Fu) = 0 — -\bnL'm{u} + kp , ( 5 . 10 ) 

which can be converted to the vector form 


kf = — XBnL' ?n{u} . 


( 5 . 11 ) 


From the results in Eqs. (5.9) and (5.11), the boundary condition matrix B appearing 
in Eq. (4.9) can be written as 


B == 


B„ O' 
L 0 0. 


( 5 . 12 ) 


For this problem, one also has that B„^ = 0. 

Having obtained expressions for the constant vector of integration, the parti- 
tioned matrices in Eqs. (4.14-4.16) are given, via Eqs. (4.10-4.12), as 


H = 
A = 
F = 


I 0 

-L'r I 

0 'fn 
.0 0 . 

(B„-I)L 

0 



( 5 . 13 ) 

( 5 . 14 ) 

( 5 . 15 ) 
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where L* can be considered as a special type of integrating matrix. (L* is actually 
the equivalent of what Vakhitov [26] refers to as a “type two” integrating matrix. 
To see the significance of this modified integrating matrix, one can write out the 
type of integration that it performs; namely, 


— r-r=/; 


(5.16) 


Other such special integrating matrices will be defined for convenience sake in 
later problems. Each one will be assigned a different subscript that will serve to 
distinguish it.) 

The matrices in Eqs. (5.13-5.15) can be used directly in Eq. (4.25) and Eq. 
(4.29) to obtain the reduced eigenvalue problem for free vibration. From Eq. (4.25), 
we obtain the result 


T = L'rLj, (5.17) 

and from Eq. (4.29) the eigenvalue problem is 

[L'rLi'fn-(l/X)I]{«} =0. (5.18) 

This result provides the eigenvalue problem for a cantilevered axial rod with variable 
stiffness and mass parameters. As can be seen, the eigenvalue problem is formed 
by specifying the integrating matrix and the corresponding discrete values for the 
diagonal matrices ' t and ' m. 

For simple problems, such as the one considered here, an alternate approach 
that often proves very convenient is to apply the integrating matrix to the form 
of the equations given in Eq. (5.1), rather than utilizing the state vector form. A 
discrete form of Eq. (5.1) can be written as 

^{Fu} + = 0. (5.19) 
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Applying the integrating matrix to Eq. (5.19) yields 

{F«} + XL'm{n} = kf . (5.20) 

Now multiplying this result through by the diagonal matrix 'r to isolate another 
derivative term gives 

+ X' 7-L' m{«} == ' rkjr . (5-21) 

Another application of the integrating matrix then produces the result 

|u} + XL' rL'm{u} = L' rk;? + ku- (5.22) 

If Eq. (5.22) is multiplied through by Bq , then use of the boundary condition «(0) = 0 
gives, as before, k„ = 0. Similarly, multiplying Eq. (5.20) by and using Fu(l) = 0 
leads to the result in Eq. (5.11). Substituting these results for the constant vectors 
into Eq. (5.22) and rearranging gives the eigenvalue problem presented in Eq. (5.18). 

From the eigenvalue problem given by Eq. (5.18), results can be obtained for 
either uniform or variable property rods. For the case of a uniform rod, the matrices 
' T and ' m are simply identity matrices. The numerical calculation of the eigenvalues 
and eigenvectors can make use of standard eigenvalue routines. For the examples 
presented, the computations were carried out on an IBM 3033. The calculations 
that form the matrix for the eigenvalue problem were performed in single precision 
arithmetic, and the eigenvalues and eigenvectors of the resulting real, nonsymmetric 
matrix were calculated with a reliable double precision routine available in the 
EISPACK eigenvalue package (see Smith [44]). Calculations were performed for both 
Jacobi and Newton integrating matrices, with the number of discretization intervals 
N varying from two to five. For each level of discretization, the maximum order 
integrating matrix was used for the given number of collocation points {N + 1). 
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The numerical solutions for axial vibration of both uniform and linearly tapered 
rods were compared with analytical solutions. For the uniform rod, analytical solu- 
tions for the constant coefficient Sturm-Liouville problem can be found in Section 
3.5 of Bisplinghoff, Ashley, and Halfman [45]. For linearly tapered rods, useful 
analytical solutions for the eigenvalues of the homogeneous differential equation are 
are given in closed form by Hildebrand and Reissner [46]. (These solutions are also 
repeated in Ref. [45], pp. 432-434.) 

In the case of the uniform rod, results for the axial frequencies obtained with 
both Newton and Jacobi integrating matrices are presented in Table 1 at the end 
of this chapter. Also presented in Table 1 are the percentage errors between the 
approximate integrating matrix results and the exact analytical solutions. Fig. 4 
gives a display of these percentage errors for the first three vibrational modes. As 
can be seen by comparing the frequency errors, both Newton and Jacobi solutions 
provide highly accurate results despite the relatively crude discretization. It is also 
clear that for those frequencies that are nearing convergence, the Jacobi solutions 
are more accurate than those provided by the Newton solutions. For the uniform 
rod, the frequencies given by the Jacobi integrating matrices always converge from 
above to the exact frequencies, unlike the oscillatory convergence of the Newton 
solutions. Since all of the available frequencies are displayed for each discretization 
level, it is useful to note how the higher frequencies behave as the discretization level 
increases. For the Jacobi solutions, it is apparent that each newly introduced modal 
frequency is quite high compared to the exact result, but rapidly converges to the 
correct solution as new discretization intervals are introduced. In comparison, each 
newly introduced Newton solution frequency is not as high, but the convergence to 
the correct result tends to be somewhat oscillatory. These convergence trends for 
uniform parameter systems are also evident in additional solutions to be presented 
in the following sections. Aside from the finer details of convergence, however, both 
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NUMBER OF INTERVALS 


Fig. 4. Error in computed frequencies, compared to exact solu- 
tions, for the first three axial modes of a uniform can- 
tilevered rod. The error for both Jacobi and Newton in- 
tegrating matrix solutions is plotted vs. the number of 
collocation intervals N. 
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Newton and Jacobi solutions have the capability to provide accurate results for the 
lower frequencies without requiring a highly refined mesh of grid points. 

For the uniform axial rod vibration just discussed, Table 2 gives the first three 
axial mode shapes. Results are shown for both Newton and Jacobi integrating 
matrix solutions based on five collocation intervals. Since the Newton solutions 
make use of evenly spaced collocation points, the Newton based mode shapes in 
Table 2 are simply presented for these collocation points without interpolation 
being required. On the other hand, the Jacobi solutions have been interpolated 
to these evenly spaced points with the same order Jacobi polynomial used in the 
solution. The results of performing such an interpolation are less than satisfactory 
since it appears that the Jacobi mode shapes in Table 2 are less accurate than 
the Newton results. But in fact, it is easily shown that the Jacobi mode shapes 
are quite accurate if evaluated at the Jacobi collocation points. This behavior 
emphasizes an important aspect of interpolation processes in general. Due to the 
fact that interpolation error is itself an oscillatory function, the highest accuracy for 
the interpolated displacement solutions can be expected to occur at the collocation 
points, and the largest error will occur between the collocation points. The same 
conclusion is reached if the Newton solutions are interpolated to the unevenly spaced 
Jacobi points. This all seems to suggest, not unexpectedly, that if mode shape 
results must be interpolated, it is probably better to use a localized least squares 
procedure based on a slightly lower order polynomial approximation. 

Vibration frequencies for a linearly tapered rod can be easily determined from 
Eq. (5.18) after introducing the following variables to describe the taper. First, the 
taper ratio of the rod is defined to be (1 — Pt), where Pt is the taper parameter. 
Each dimension in the cross section of the rod varies linearly as (1 — pix). For this 
description of the taper, it follows that the mass and stiffness parameters m and 
EA appearing in Eq. (5.2) can be written as 
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Fig. 5. Error in computed frequencies, compared to exact solu- 
tions, for the first three axial modes of a linearly tapered 
cantilevered rod with a taper ratio of one half (/?( = 0.5). 
The error for both Jacobi and Newton integrating matrix 
solutions is plotted vs. the number of collocation inter- 
vals N. 
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m = mjt(l - 


(5.23) 


EA=={EAUl-0tx)\ 

Using these definitions for the mass and stiffness, the frequencies of a linearly 
tapered rod with = 0-5 are given in Table 3. The percentage errors between 
the exact solutions and the integrating matrix solutions are also displayed in Fig. 5 
for the first three vibrational modes. As evident from the eigenvalue error results, 
neither the Jacobi nor the Newton solutions converge quite as rapidly as in the 
uniform rod case, which is to be expected. For this linearly tapered rod, the Jacobi 
results do not necessarily converge to the exact frequencies from above as they do for 
uniform parameter problems. Nonetheless, the convergence of the Jacobi solutions 
is quite predictable in that all modes, except the first, appear at a high frequency 
and only very slightly overshoot on the negative side enroute to converging. The 
Jacobi solutions still provide faster convergence to the lower modes. In contrast, the 
Newton solutions retain their oscillatory convergence character, with the frequency 
error changing sign for each added collocation interval. Aside from this oscillatory 
behavior, the magnitude of the Newton solution errors can be considered quite 
small. It is well known that one cannot always guarantee convergence from above 
for eigenvalues obtained from a collocation procedure. The fact that there are sign 
changes in the error of both Jacobi and Newton frequencies for this nonuniform 
parameter problem is indicative of the collocation nature of the integrating matrix 
solutions. 

5.1.2 Rods with Discontinuous Mass and Stiffness 

Because of the common occurrence of structural elements with parameter 
discontinuities, it is useful to examine the hybrid state vector technique as applied to 
axial vibration of rods having stepwise jumps in mass and stiffness. The eigenvalue 
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problem for a discontinuous rod is still given by Eq. (5.18), which was originally 
derived for continuous parameter rods. The only difference arises in the proper 
specification of the integrating matrix and the diagonal coefficient matrices 'r and 
'm. The methods for specifying the integrating matrix and the corresponding 
function values for a discontinuous integrand were discussed in Section 3.2. 

To allow comparison with alternate solutions, the problem chosen for discon- 
tinuous axial rod vibration is the same as a problem discussed by Hodges [47], in 
which solutions were obtained by a Ritz method. The cantilevered rod, which is 
normalized to unit length and clamped at tt = 0, is divided into three segments. 
The first segment extends from x = 0 to * = 0.25, the second segment from 0.25 
to 0.75, and the third segment from 0.75 to 1.0. Mass and stiffness parameters in 
the first and third segments are assigned the same values, and within each segment 
the parameters are assumed constant. The ratios between the parameters of the 
different segments are given by 


{EA)2 

{EA)i T2 ’ mi 


(5.24) 


where 




1+7 


_ _ 2 
mi = mz — — -3 

1 + e 


(5.25) 


Thus, when 7 = 1 and $ = I there is no discontinuity in the parameters. 

Integrating matrix solutions are obtained for the discontinuous case by cal- 
culating an integrating matrix for the piecewise integration of the three segments. 
Maximum order weighting matrices are assumed in this calculation. An equal num- 
ber of collocation intervals are used for each segment; thus, it is appropriate to 
describe the solution in terms of the number of intervals per segment. 
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For 7 = 10 and ^ = 100 the results are presented in Table 4 for the axial 
frequencies obtained with both Jacobi and Newton integrating matrices. Fig. 6 
displays for the first three vibrational modes the percentage error between the 
exact solutions (Ref. [47]) and the integrating matrix solutions. Similar to the 
earlier discovery for uniform axial rod vibration, the Jacobi solutions for a piecewise 
uniform rod offer the advantage of a rapid, monotonic convergence to the exact 
result from above. From a practical viewpoint, either the Jacobi or the Newton 
integrating matrix yields highly accurate results without requiring an extremely 
fine mesh of grid points. 

An interesting point to notice in Fig. 6 is that the eigenvalue error for the 
fundamental mode bottoms out and then begins to rise slightly. The reason for this 
occurrence is that the numerical precision limit has been reached for the calculation 
of this eigenvalue. As a consequence, the last digits of the result change rather than 
remaining exactly at the converged value. D.H. Hodges (in a private communication) 
stated that this same phenomenon was encountered when obtaining similar high 
accuracy solutions with a variable order Ritz finite element method. The apparent 
cure for this is to increase the precision level of the computations. 

Because integrating matrices with end points were used to solve this problem, 
one should note that two identically located collocation points are present at each 
discontinuity. Since the solution values will be identical for each of these points, 
a constraint is available that can be used to eliminate the solution degrees of 
freedom at one of the points. The rectangular transformation matrix resulting 
from this constraint can be applied as discussed in Appendix F. It turns out that 
this transformation matrix, which consists only of ones and zeroes, produces a 
diagonal matrix when premultiplied by its transpose. Therefore, the pseudoinverse 
mentioned in Appendix F can be obtained quite simply. It should be pointed out, 
however, that for small sized problems it is not necessary, and generally not worth 
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INTERVALS PER SEGMENT 


Fig. 6. Error in computed frequencies, compared to exact solu- 
tions, for the first three axial modes of a cantilevered rod 
with discontinuous stiffness and mass (7 = 10 ; 0 = 100 ). 
The rod has three uniform segments between the points 
= 0 , .25, .75, 1 . The error for the integrating matrix 
solutions is plotted vs. the number of collocation inter- 
vals per segment N. 
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the trouble, to apply this constraint. The exception would be a large sized problem. 
In that case one might expect a significant reduction in the number of degrees of 
freedom present in the solution. 

A final note concerning axial vibration of a discontinuous rod is that solutions 
were also tried in which the discontinuities were ignored. That is, an integrat- 
ing matrix was used which allowed the integrations (i.e., the interpolations) to 
proceed across the discontinuity boundaries. As expected, the eigenvalue solutions 
approached the exact results when enough collocation points were taken. Although 
a reasonably accurate solution for the fundamental mode could be obtained without 
too much effort, matching the accuracies of Table 4 for the higher modes required a 
somewhat larger number of collocation points. Furthermore, the convergence with 
this approach was highly oscillatory for some modes and seemed unpredictable for 
anything but the fundamental mode. 


5.1.3 Rods with Elastic Restraint 

A type of boundary condition that occurs quite often is that of elastic restraint. 
Methods for handling this type of boundary condition are easily demonstrated for 
a uniform axial rod. It will be assumed that the axial modes of vibration are to be 
found for a rod that is cantilevered at one end and has an axial spring restraint at the 
opposite end. Fig. 7 illustrates a rod with these boundary conditions. Since the axial 
restraining spring applies a force to the rod that is proportional to the displacement 
at the end of the rod, the boundary condition can be included as a displacement 
dependent loading in the Kfd submatrix of Eq. (4.15). In the present example, 
however, the restraint will be introduced instead through the nonhomogeneous 
boundary condition matrix 

The nondimensional equations describing the rod are given by Eqs. (5. 3-5. 4) 
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Fig. 7. A cantilevered axial rod with a spring restraint boundary 
condition 

and Eqs. (5. 6-5. 8). Boundary conditions for the rod are 


«(0) = 0, F„(l) = -Argu(l), 


( 5 . 26 ) 


where 


Hg — 


kg 

{EA)k • 


( 5 . 27 ) 


As the nondimensional restraining spring stiffness kg varies from zero to infinity, 
the boundary conditions for the rod vary from cantilevered to clamped-clamped. 

Solutions for this problem are obtained by the same approach as presented in 
Section 5.1.1. In fact, the only changes to be made are to the boundary condition 
solution in Eqs. (5.10-5.11) and to the H matrix in Eq. (5.13). The boundary 
condition solution for the spring restrained end is given by 


bn{Eu} = -Tcgbn{u} = ~Xb^L'm{«) + kp , ( 5 . 28 ) 


which yields the vector form 

kf = XBnL'rn{«} — FBBn{u) . ( 5 . 29 ) 
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From this result, one finds that the boundary condition matrix B is still given by 
Eq. (5.12), and the nonhomogeneous boundary condition matrix is now given by 


Bfi/i = 


0 ^gBfi 

LO 0 . 


(5.30) 


With the definition for H in Eq. (4.11), the replacement for Eq. (5.13) is found to 
be 


H = 


I 


Tcb Bfi 


(5.31) 


L-L'r I J 

If one now makes use of Eq. (4.27), along with Eqs. (5.14^5.15), the eigenvalue 
problem can be written as 


[I*L'rLt'rn-(l/X)I]{u} =0 


(5.32) 


where 


I* = [I + A:gL'rB„]-^ 


(5.33) 


The eigenvalue problem in Eq. (5.32) applies to vibration of a nonuniform rod and 
can be solved for any value of restraining spring stiffness. 

Although it might appear to be a formidable task, it is actually quite easy to 
analytically invert the expression in Eq. (5.33), and thus obtain a simpler form for 
the eigenvalue problem. The reason that Eq. (5.33) can be easily simplified is that 
B„ contains unit terms only in the last row, and zeroes elsewhere. Thus, for all 
practical purposes, the second expression in the summation of Eq. (5.33) yields a 
vector. This allows easy inversion of the entire matrix when applying inversion by 
matrix partitioning. In general terms, one can write 


[I + FgL'rB„] 


•I ei 
.0 C2. 


(5.34) 


- 77 - 


where the matrix has been partitioned such that C 2 is a scalar and ei is a vector. 
The inverse of Eq. (5.34) is 


I* 


'I -ei/c2' 

.0 lfc2 . 


This result can be conveniently written in the form 


I* = [l-vksVrBn] 


( 5 . 35 ) 


( 5 . 36 ) 


where 


_ _ 1 

^2 i + k,bj;vri' 


( 5 . 37 ) 


For the case of a uniform rod (' r = 'm = I), the eigenvalue problem given by 
Eq. (5.32), in conjunction with Eq. (5.36), can be written as 


[rLLt-(l/X)I]{«}=0 


( 5 . 38 ) 


where 




I- 


1 + if, 


LB„ 


( 5 . 39 ) 


Note that in obtaining Eq. (5.39) from Eq. (5.37) a special identity was used; namely, 


biLi = 1. 


( 5 . 40 ) 


This identity applies when the integrating matrix has been written for a normalized 
interval of unit length. It is well worth remembering this identity since it appears 
often in solving for various types of boundary conditions. 

Exact eigenvalues for the uniform cantilevered rod with spring restraint are 
easily found by analytically solving the differential equations. It can be shown that 
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the eigenvalues are given by the roots of the transcendental equation 

Jfg tan i/X + \/X = 0, (5.41) 

and the corresponding axial mode shapes are obtained from 

u{z) = C sin\/\x. (5.42) 

For the exact solutions, as Fg-»0, \/X— with n = 1,3, 5, . . . , oo. This is the 
result for a purely cantilevered rod. On the other hand, as ►oo, \/X— ►nn-, with 
n = 0, 1 , 2, . . . , oo, which is the solution for a rod clamped at both ends. 

Presented in Table 5 for a uniform rod are nondimensional frequencies calcu- 
lated from Eq. (5.38). Comparisons are made with the exact results of Eq. (5.41). 
The frequencies are presented for the first three modes and are calculated for several 
values of the restraint spring stiffness kg. These calculations make use of a Jacobi 
integrating matrix corresponding to six collocation points. The exact frequencies 
listed under = 1 X 10® actually correspond to kg = oo. This gives a valid com- 
parison with the integrating matrix solutions since values of the spring stiffness 
larger than 1 X 10® produce no significant changes in the integrating matrix results 
when using single precision calculations. 

5.1.4 Rods with Concentrated Mass 

Another important example deals with the vibration of rods having concen- 
trated mass points in addition to continuously distributed mass properties. The 
approach to solving problems with concentrated masses involves a simple extension 
of methods already presented for continuous parameter problems. This extension 
of the foregoing methods is presented in Appendix E, which discusses the treatment 
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of any type of concentrated loading. For the problem at hand, the point loads arise 
from concentrated masses. 

The problem to be considered in this section is the axial vibration of a can- 
tilevered rod with a tip mass. The method of solution will be that of Appendix 
E, but one should note that for this particular problem the mass provides a tip 
boundary condition that can also be handled with the approach given in the pre- 
vious section. Since the solution turns out to be an extension of that presented in 
Section 5.1.1, the matrices defined in Eqs. (5.13-5.15) still apply. To account for 
the concentrated inertia loads, however, it is necessary to include from Eq. (E.ll) 
the matrix 


1 

o 

7 

e 


rsf 01 

1 0 -d 


1 

o 


(5.43) 


Notice that this matrix is the same as F in Eq. (5.15) with the summing matrix S 
substituted for the integrating matrix. A second matrix that must be added is the 
matrix containing the concentrated mass terms specified at the collocation points. 
For the problem under consideration this matrix is 


A+ 




m 




(5.44) 


where = mnlniRl is the tip mass nondimensionalized by the total mass of 
a reference rod (i.e., is a total mass and rriR is a running mass). With the 
matrices defined in Eqs. (5.43-5.44), the eigenvalue problem given by Eq. (E.16) can 
be written as 


[(L'rLt'm-hL'rSrm;t)-(lA)*]{«} =0. (5.45) 
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This eigenvalue problem is the extension of Eq. (5.18) to include point mass terms. 
For a uniform beam, Eq. (5.45) simplifies to 

[(LLt + LSt'm+)-(l/X)I]{ti} = 0 (5.46) 

where it should be noted that the lumped tip mass will only affect terms in the last 
column of the matrix from which the eigenvalues are found. 

Exact solutions for the eigenvalues of a similar Sturm-Liouville problem can 
be found in Section 5.9 of Meirovitch [42]. The transcendental equation for the 
eigenvalues is given by 


\/X tan — 1/m^ = 0. (5.47) 

Upon examining this equation, one finds that when mn— *^0, then \/X— »-n7r/2, with 
n = 1, 3, 5, . . . , oo. Similarly, as m„ — ► oo, y/\ — ► nn, with n = 0, 1, 2, . . . , oo. In physical 
terms, as the tip mass tends to zero, the eigenvalue solution reduces to that for a 
simple cantilever rod. As the tip mass grows very large, the solution approaches 
that for a rod clamped at both ends. 

Table 6 gives for a uniform rod the nondimensional frequencies calculated 
from Eq. (5.46). Comparisons are made with the exact results of Eq. (5.47). The 
frequencies are given for the first three modes and are calculated for several values 
of the tip mass ratio mj}". Calculations are based on a Jacobi integrating matrix 
corresponding to six collocation points. The exact frequencies that are listed under 
m'^ = 1 X 10^ actually correspond to = oo. This gives a valid comparison 
with the integrating matrix solutions since values of the tip mass ratio larger than 
1 X 10* produce no further changes in the integrating matrix results when using 
single precision calculations. 
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5.2 Bending Vibration of Beams 

Applications of the hybrid state vector approach are now examined for beam 
vibration problems. Since beam vibration analyses are so common in structural 
dynamics and aeroelasticity, it is useful to give them a separate treatment. The 
following problems demonstrate the application of hybrid state vector techniques to 
the lateral vibration of beams with various type of boundary conditions, including 
unrestrained beams possessing rigid body modes. Sample results are presented for 
uniform beams. 

The differential equations that describe lateral vibration of isotropic beams are 
well known. Rather than rederiving the nondimension al state vector equations, 
they can be extracted from the anisotropic beam equations given in Section 2.3. 
After including inertia terms and neglecting transverse shear effects, the appropriate 
equations can be written as 



M 


• 0 

1 

0 

O' 


M' 


'0 

0 

0 

O' 


'M' 

d 

V 


0 

0 

0 

0 


V 


0 

0 

0 

m 


V 

— { 


> = 







i- X 







dx 

1 


m 

0 

0 

0 


7 


0 

0 

0 

0 


7 


‘ w - 


. 0 

0 

-1 

0. 


* W ' 


.0 

0 

0 

0. 


. w ‘ 


where 


EI = 


{EI)r 
E l ’ 


m = 


m 

mu ’ 


and 


{EI)r 


(5.48) 


(5.49) 


(5.50) 


Other nondimensional terms are as defined in Eq. (2.29) and Eq. (2.31), with c taken 
to have unit value for a beam. 
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Irrespective of boundary conditions, one can begin the solution of Eq. (5.48) by 
discretizing, multiplying through by the global integrating matrix, and adding the 
constant vector of integration. In a manner analogous to Eqs. (5. 6-5. 8 ), this result 
can be expanded into the set of equations 


{M} = L{V} -Mcn, 

(5.51) 

{F} = -h k„ 

(5.52) 

{ 7 } =L'M{M}+k^ 

(5.53) 

{«)} = -L{ 7 } + k„, . 

(5.54) 


Boundary conditions can be applied to these equations to solve for the constant 
vectors of integration. 

The basic approach to be used in solving for the constants of integration of 
two-point boundary value problems is demonstrated in Section 5.1.1. Only the 
important aspects of boundary condition determination will be addressed in detail 
in the following problems. It is worth noting, however, that there is a recommended 
procedure to simplify the solution for boundary conditions. One should first apply 
any homogeneous or nonhomogeneous boundary conditions at the end * — 0 . It is 
assumed that this end is the starting point for the integrations performed by the 
integrating matrix. The second step is to apply boundary conditions at the point 
* = 1 (assuming, of course, a normalized interval * = [ 0 , 1 ] ). This second step often 
involves applying boundary conditions to state variables that were already used 
in the boundary conditions at ar = 0 . When this occurs, one must express these 
variables in terms of the remaining unused variables to insure that all unknown 
constants are determined. An important requirement when solving for constants of 
integration is that the resulting boundary condition matrix B be block diagonal. 
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5.2.1 Cantilevered Beam 


A cantilevered beam offers simple boundary conditions from which the con- 
stants of integration are readily determined from Eqs. (5.51-5.54). The boundary 
conditions to be satisfied are given by 

1(0) = 0, W(0) = 0, M(1) = 0, F(1) = 0. (5.55) 

It is clear from the first two conditions that = 0. The boundary conditions 

at the free end of the beam yield 


km = -B„L{F} 
kv — XBnL'm{w} . 

With these constants determined, one finds that 


(5.56) 

(5.57) 


1 

c 

CD 

0 

0 

0 ■ 


r 

0 

0 

0 ■ 

0 

(Bn-I)L 

0 

0 


0 

Lf 

0 

0 

0 

0 

-L 

0 


0 

0 

-L 

0 

0 

0 

0 

-L. 


.0 

0 

0 

-L. 


(5.58) 


and as a result, 


■ I Lf 0 O' 

0 10 0 

H = 

-L':^ 0 I 0 

.0 0 L I. 

From Eq. (4.25), one then obtains 

' L'mLl[ -L'mLf 

T = 


(5.59) 


(5.60) 
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which, when substituted into Eq. (4.29), yields the eigenvalue problem 

[L2's7Ll^'m-(l/X)I]{w} = 0. (5.61) 

This simplifies for a uniform beam since 'IT and 'm are then identity matrices. 
It should be noted that the row and column of Eq. (5.61) corresponding to the 
constrained lateral displacement at the fixed end can be deleted before solving for 
the eigenvalues and eigenvectors. 

Results are presented in Table 7 for the nondimensional frequencies of a uniform 
cantilevered beam. The numerically computed results are based on a Jacobi in- 
tegrating matrix for five collocation intervals and are compared with exact results 
tabulated by Hunter [27], One can also refer to Hunter for detailed Newton in- 
tegrating matrix solutions of the cantilevered beam problem. 

Mode shapes corresponding to the first three modes of the Jacobi solutions 
are given in Table 8. The accuracy of the higher modes can be improved further 
by adding rotary inertia and transverse shear terms to the analysis. The mode 
shapes in Table 8 were interpolated to the evenly spaced grid points with a Jacobi 
polynomial of the same order used in the integrating matrix solution. Some of the 
aspects of this type of interpolation are discussed in Section 5.1.1. 

5.2.2 Simply Supported Beam 

A simply supported beam requires boundary conditions to be applied to both 
moment and lateral displacement variables at both ends of the beam. The familiar 
boundary conditions for a simply supported beam are concisely stated as 

M(0) = 0, tu(0) = 0, M(1) = 0, W(1)=0. (5.62) 

The boundary conditions atx = 0 applied toEqs. (5.51) and (5.54) yield immediately 

SjTj = = 0 . (5.63) 
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To apply boundary conditions at x = 1, it is necessary to write expressions for M 
and w in terms of the other unknown constants appearing in Eqs. (5.52) and (5.53). 
These expressions are 


{M} = L(— XL'm){w} + Lki, (5.64) 

{w} = -L(L':^){M}«; - Lk.y. (5.65) 


After applying the boundary conditions at s = 1 and making use of the two 
identities ky = Ikj and b^Ll = 1, one obtains the results 


kt, = — BnL(XL' m){w} 

k^ = -B„L(L'E7){M), 


which leads to 


(5.66) 


•-L 

0 

0 

0 • 


■-L 

0 

0 

0 ■ 

0 

(B„L-I)L 

0 

0 


0 

LI 

0 

0 

0 

0 

(B„L-I)L 

0 


0 

0 

L| 

0 

. 0 

0 

0 

-L. 


. 0 

0 

0 

-L. 


(5.67) 


From Eq. (5.67), one then obtains in the manner of Section 5.2.1 that 


T = 


i-l'EIL 
"2 


L-LLr£?/L 


V.Ll'E~ILLl. 


from which follows the reduced eigenvalue problem 


(5.68) 


(LL^'E/LL|'m-(l/X)I]{w) =0. 


(5.69) 


Results are presented in Table 9 for the non dimensional frequencies of a uniform 
simply supported beam. The numerically computed results are based on Jacobi 
integrating matrices and are compared with exact results given in Ref. [45]. 
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5.2.3 Hinged-Free Beam 


The hinged-free beam considered in the following analysis is assumed to be 
simply supported at x = 0 and unrestrained at af = 1; this leads to a single rigid 
body mode in rotation. Boundary conditions for this beam are specified by 


M(0) = 0, W(0) = 0, M(1) = 0, y(l) = 0. (5.70) 


From the conditions at -p = 0 it is easily found that 

kjTi = k«; = 0, (5.71) 


and at 2 = 1 one readily has 


k„ == XB„L' m{w} . (5.72) 

To complete the solution for the remaining constant of integration, it is necessary 
to write an expression for M in terms of k-y. After substituting Eqs. (5.71-5.72) into 
Eqs. (5.51-5.54), this expression can be obtained in the form 

(M) = -XLLi'mL(L'E7{M} k^). (5.73) 

From the boundary condition M(l) = 0, one then obtains the result 

k.y = -c.yB„LLt'mL(L'E7){70'} (5.74) 


where 


bJ^LLj'mLl 


(5.75) 


- 87 - 



The foregoing solutions for the constants of integration lead to the expressions 


•-L 

0 
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0 ■ 


■-L 
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(c^B„LLt'mL-I)L 

0 


0 

0 


0 

. 0 

0 

0 
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-L. 


(5.76) 


and 


T = 


■ Ll'EIL 
-LL|'m 


-LfEILLl 

LLfmLL*i. 


(5.77) 


From Eq. (5.77) follows the reduced eigenvalue problem 


[LL3'F;/LL*'m-(l/X)I]{w} =0. 


(5.78) 


Results are presented in Table 10 for the nondimensional frequencies of a 
uniform hinged-free beam. The numerically computed results are based on Jacobi 
integrating matrices and are compared with exact results obtained from the expres- 
sion tan - tanh 0 = 0, where \/\ = 0^ is the nondimensional frequency. 

Solutions for the eigenvalues of Eq. (5.78) can be obtained without having to 
apply constraint equations. Constraint equations do exist, however, in problems 
with rigid body modes. The presence of constraints in an eigenvalue problem will be 
evident from the existence of zero eigenvalues. With finite precision calculations, the 
zero eigenvalues are not guaranteed to be exactly zero for nonsymmetric eigenvalue 
problems like that in Eq. (5.78). In practice, these zero eigenvalues can appear as 
very small positive or negative quantities, and in some cases, complex quantities 
with very small real parts. The appearance of these parasitic eigenvalues is not 
an indication that a solution has gone astray; parasitic eigenvalues can simply be 
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neglected. In Eq. (5.78) the eigenvalues are related to reciprocals of the vibration 
frequency; therefore, the parasitic eigenvalues will equivalently manifest themselves 
as very high frequencies. 

For problems with rigid body modes, conservation of linear and angular momen- 
tum provide the constraints that can eliminate associated zero eigenvalues. The 
constraint arising from linear momentum conservation is 

wmdx = 0, (5.79) 

which, for free vibration, leads to the discrete constraint equation 

b^L{tf;rn} = 0. (5.80) 

Similarly, the constraint given by angular momentum conservation is 

ibmxdx = 0, (5.81) 




which results in 


bJ[L{«;mar} = 0. (5.82) 

These constraint equations can be used to reduce an eigenvalue problem as discussed 
in Appendix F. 


5.2.4 Free-Free Beam 

A free-free beam possessing rigid body translation and rotation modes presents 
a slightly more complicated solution for the constants of integration. Nonetheless, 
the procedure to be followed is similar to that used in previous problems. The 
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boundary conditions for the free-free beam are given by 


M(0) = 0, y(0) = 0, M(1) = 0, y(l) = 0. (5.83) 


The application of these boundary conditions to Eqs. (5.51-5.54) gives the constants 
of integration 


— 0 
=0 

kty = rn(— L){7} 

= -c^B„L2'7nL|(L'M){M} , 


(5.84) 


where 


Cw; 


b^L'ml 


(5.85) 


C/v 


b^:L2'mL|i 


and 


L| =(c^B„L'm-I)L. 


(5.86) 


These results then lead to the expressions 
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L4J 


(5.87) 


and 

r L|'E7L L|'E7 l2 ■ 
L|L|'E7L L|L|'E7l2 


(5.88) 
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With the aid of Eq. (5.88), the reduced eigenvalue problem is given by 


{qLrm2'm-(i/x)i]{w} =0. 


For a uniform property beam the eigenvalue problem reduces to 


[l|l|l2-(iA)i]{w}-o 


where 


(5.89) 


(5.90) 


Ll = (c^B„l2l| - I)L 


(5.90) 


and 


b^L2L|i 


(5.91) 


Note that is the same as given in Section 5.2.2. 

Table 11 gives the results for the nondimensional frequencies of a uniform 
free-free beam. The numerically computed results are based on Jacobi integrating 
matrices and are compared with exact results obtained from (cosh ,3) cos ^ — 1=0, 
where is the nondimensional frequency. If desired, constraint equations 

that reduce out the rigid body modes can be applied as discussed in the previous 
section. 


5.3 Buckling of a Rotating Beam 

In the following analysis, the state vector approach will be used to obtain 
solutions for buckling instabilities of the cantilevered, inward oriented rotating beam 
pictured in Fig. 8. Buckling is assumed to take place in the plane of rotation. 
Since this problem has been thoroughly examined by many investigators, results 
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Fig. 8. Rotating beam geometry 

are available with which to compare the integrating matrix solutions. This analysis 
demonstrates that solutions can be easily obtained for many problems involving 
aeroelasticity, buckling, or vibration of rotating machinery. 

The state vector equations to be employed in the buckling analysis compare 
directly with the equations used by Steele and Barry [13]. Considering bending 
deformations in the plane of rotation, the nondimensional form of these equations 
is given by 
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(5.92) 


where the nondimensional tension parameter Ft is given as 

f,(x) = (l/2)|(c-l)2-(a,-l)2l 

with 


(5.93) 


tts = 


(5.94) 
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I 


The nondimensional spin parameter Xg is given by 


‘ {ei)r • 


(5.95) 


Except for the displacement dependent load terms, these equations are very similar 
to those presented in Eq. (5.48). The solution to Eq. (5.92) is readily obtained by 
using the T matrix presented in Eq. (5.60) for cantilever boundary conditions. The 
resulting buckling eigenvalue problem is given by 


where 


[TA;.^-(l/Xg)I] 



= 0 


T^fd = 


-L'EIL’l'Ft 


-L'EILf'm' 

L^'EILf'm_ 


(5.96) 


(5.97) 


Results are given in Table 12 for the nondimensional critical buckling speed 12 
of a uniform beam, where 12 = af\A7- The buckling speeds, which are presented 
for varying values of the parameter ag, are compared with numerical solutions 
given in Ref. [13]. These comparison solutions consist of results from an asymptotic 
iteration scheme as well as results from a high order Ritz finite element method 
due to Hodges [48]. Further elaboration on these critical buckling speed solutions 
is given by White, Kvaternik, and Kaza [31] and by Peters and Hodges [49]. 


5.4 Deflection of Beams 

Static beam deflections and forces can be easily and accurately calculated 
with the hybrid state vector approach. The nondimensional state vector equations 
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describing lateral deflection of a beam under an applied loading p are 


where 
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(5.98) 


Wh' 


(5.99) 


Solutions to Eq. (5.98) are given by Eqs. (4.23) and (4.24). If aerodynamic loads 
are not present, one must take X — 0. Noting that = 0, then for cantilever 
boundary conditions the deflection and force solutions are given by 


yz) = (5.100) 

yf = H^^L|arjp* (5.101) 

By using the T matrix in Eq. (5.60), these solutions can be further refined to give 


and 



(5.102) 


(5.103) 


As an indicator of solution accuracy, one can easily verify with a sample calcula- 
tion that a four interval Jacobi integrating matrix solution will provide the exact 
deflections for a uniform cantilevered beam having a constant load distribution. 


- 94 - 



It should be noted that matrix inversions are not required in solving for either 
forces or displacements. Deflection solutions similar to these are used to obtain 
static aeroelastic lift distributions in the next chapter. 
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Table 1 


Nondimensional frequencies, for the axial modes of a uniform 
cantilevered rod with N collocation intervals. (X = I{EA)r) 


N 

Mode 

Jacobi 

Newton 

Exact 

% Error 

Jacobi Newton 

2 

1 

1.582576 

1.582576 

1.570796 

.749938 

.749938 


2 

7.582580 

7.582580 

4.712389 

60.9073 

60.9073 

3 

1 

1.571009 

1.569703 

1.570796 

.013560 

-.069583 


2 

4.962687 

4.826099 

4.712389 

5.31149 

2.41300 


3 

15.39166 

14.25626 

7.853982 

95.9727 

81.5163 

4 

1 

1.570799 

1.570732 

1.570796 

.000191 

-.004074 


2 

4.735692 

4.665351 

4.712389 

.494505 

-.998177 


3 

8.807365 

8.070704 

7.853982 

12.1388 

2.75939 


4 

25.64255 

21.64282 

10.99557 

133.208 

96.8321 

5 

1 

1.570797 

1.570803 

1.570796 

.000064 

.000446 


2 

4.713869 

4.706141 

4.712389 

.031407 

-.132587 


3 

8.011505 

7.631141 

7.853982 

2.00564 

-2.83730 


4 

13.27345 

11.23810 

10.99557 

20.7163 

2.20567 


5 

38.40523 

29.57599 

14.13717 

171.661 

109.207 


Table 2 


Axial vibration mode shapes, {«}, for a uniform cantilevered rod 
(Numerical solutions obtained with five collocation intervals) 


X 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

Mode 1 

Jacobi 

0.00000 

0.30903 

0.58779 

0.80902 

0.95107 

1.00000 

Newton 

0.00000 

0.30902 

0.58779 

0.80902 

0.95106 

1.00000 

Exact 

0.00000 

0.30902 

0.58779 

0.80902 

0.95106 

1.00000 

Mode 2 

Jacobi 

0.00000 

-.81474 

-.94453 

-.30816 

0.58191 

1.00000 

Newton 

0.00000 

-.81516 

-.95512 

-.31182 

0.59050 

1.00000 

Exact 

0.00000 

-.80902 

-.95106 

-.30902 

0.58779 

1.00000 

Mode 3 

Jacobi 

0.00000 

0.91354 

-.01155 

-.87489 

-.09662 

1.00000 

Newton 

0.00000 

0.98697 

-.03503 

-.95027 

-.08023 

1.00000 

Exact 

0.00000 

1.00000 

0.00000 

-1.0000 

0.00000 

1.00000 
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Table 3 


Nondimensional frequencies, for the axial modes of a linearly 
tapered cantilevered rod with N collocation intervals. [0t = 0.5 ; 
X = I{EA)r] 


N 

Mode 

Jacobi 

Newton 

Exact 

% Error 

Jacobi Newton 

2 

1 

1.613095 

1.813095 

1.652805 

-2.40258 

-2.40258 


2 

11.15869 

11.15869 

3.627852 

207.584 

207.584 

3 

1 

1.643238 

1.665064 

1.652805 

-.578834 

.741709 


2 

3.587127 

3.535655 

3.627852 

-1.12257 

-2.54137 


3 

22.39381 

20.38371 

5.807501 

285.601 

250.989 

4 

1 

1.652913 

1.649068 

1.652805 

.006534 

-.226100 


2 

3.580021 

3.658788 

3.627852 

-1.31844 

.852736 


3 

6.040474 

5.622681 

5.807501 

4.01159 

-3.18278 


4 

37.76820 

30.95097 

8.215935 

359.694 

276.719 

5 

■■ 

1.652826 

1.653185 

1.652805 

.001271 

.022991 


WM 

3.633674 

3.608995 

3.627852 

.160481 

-.519784 


H 

5.721657 

5.887524 

5.807501 

-1.47816 

1.37793 


B 

9.001691 

7.591113 

8.215935 

9.56381 

-7.60500 


mm 

57.31683 

42.52789 

10.27790 

457.671 

313.780 
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Table 4 


Nondimensional frequencies, for the axial modes of a three 
segment cantilevered rod with discontinuous stiffness and mass. 
There are N collocation intervals per segment. (X = 

=7=10,?= 100; Error based on sLx significant figures) 


N 

Mode 

Jacobi 

Newton 

Exact 

Cfr 

/o 

Jacobi 

Error 

Newton 

■I 

1 

.8336431 

.8336431 

.8267091 

.838747 

.838747 


2 

22.20007 

22.20007 

6.097328 

264.095 

264.095 

■ 

3 

121.7023 

121.7023 

11.98907 

915.105 

915.105 

B 

1 

.8267307 

.8267298 

.8267091 

.002661 

.002540 

■ 

2 

6.729430 

6.729422 

6.097328 

10.3668 

10.3667 


3 

19.03003 

19.03001 

11.98907 

58.7275 

58.7275 


4 

41.08241 

41.08287 

17.57944 

133.696 

133.699 

■ 

5 

85.56076 

85.56099 

19.46456 

339.571 

339.572 

3 

1 

.8267099 

.8267071 

.8267091 

.000121 

-.000242 


2 

6.139068 

6.082418 

6.097328 

.684562 

-.244533 


3 

14.52343 

13.85007 

11.98907 

21.1384 

15.5224 


4 

19.45551 

19.33504 

17.57944 

10.6722 

9.98609 


5 

38.17844 

37.79570 

19.46456 

96.1427 

94.1766 

4 

1 

.8267105 

.8267104 

.8267091 

.000121 

.000121 


2 

6.099150 

6.085790 

6.097328 

.029849 

-.189263 


3 

12.34547 

11.82625 

11.98907 

2.97270 

-1.35790 


4 

18.83118 

18.54741 

17.57944 

7.12083 

5.50645 


5 

25.28983 

22.48436 

19.46456 

29.9271 

15.5143 

5 

1 

.8267107 

.8267107 

.8267091 

.000242 

.000242 


2 

6.097375 

6.097551 

6.097328 

.000820 

.003608 


3 

12.03000 

11.86734 

11.98907 

.341143 

-1.01592 


4 

18.17894 

17.10370 

17.57944 

3.41024 

-2.70601 


5 

20.28574 

19.38240 

19.46456 

4.21843 

-.422331 
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Table 5 


Nondimensional frequencies, for the axial modes of a uniform 
cantilevered rod with variable stiffness elastic restraint at the free 
end. Results were obtained with a Jacobi integrating matrix using 
six collocation points. = k}[EA)R ; X = J^mR^j[EA)R) 


Mode 


He = 0.25 


1 

(exact) 

1.57080 

(1.57080) 

1.71551 

(1.71551) 

3.14161 

(3.14159) 

2 

(exact) 

4.71387 

(4.71239) 

4.76044 

(4.76481) 

6.30594 

(6.28318) 

(exact) 

8.01150 

(7.85398) 

8.04779 

(7.88567) 

10.1059 

(9.42478) 


Table 6 

Nondimensional frequencies, for the axial modes of a uniform 
cantilevered rod with tip mass. Results were obtained with a Jacobi 
integrating matrix using six collocation points, (mjj' = mnjmRt ; 
X = u?mRej{EA)R) 


Mode 



mjj- = 1.0E8 


1.57080 

(1.57080) 

1.20770 

(1.26459) 

0.00010 

(0.00000) 


4.71387 

(4.71239) 

3.98864 

(3.93517) 

3.24838 

(3.14159) 

3 

(exact) 

8.01150 

(7.85398) 

7.00118 

(6.81401) 

6.52342 

(6.28319) 




























Table 7 


NondimensioDa.1 bending frequencies, of a uniform cantilever 
beam for a Jacobi integrating matrix solution using five collocation 
intervals. (\ = j{EI)ii) 


Mode 

. ' ^ 

2 

3 

4 

5 

Calculated 

3.516022 

22.04884 

64.15257 

178.0921 

I032.3I9 

E.xact 

3.516015 j 

22.03449 

61.69721 

120.9019 

199.8595 


Table 8 


Lateral bending mode shapes, {«;}, of a uniform cantilever beam 
(Jacobi integrating matrix solution using five collocation intervals) 


Mode\3 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

1 

(exact) 

0.00000 

(0.00000) 

0.06386 

(0.06387) 

0.22990 

(0.22989) 

0.46114 

(0.46114) 

0.72547 

(0.72548) 

1.00000 

(1.00000) 

2 

(exact) 

0.00000 

(0.00000) 

-.30660 

(-.30106) 

-.68022 

(-.68347) 

-.58537 

(-.58948) 

0.06418 

(0.07004) 

rTooooo 

(1.00000) 

3 

(exact) 

0.00000 

(0.00000) 

0.60087 

(0.60450) 

0.45954 

(0.52593) 

-.40455 

(-.47377) 

-.41063 

(-.39488) 

1.00000 

(1.00000) 
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Table 9 


Nondimensional bending frequencies, of a uniform, simply 

supported beam with N collocation intervals; (X = w ^mRt^l{EI)R) 


N 

Mode 

Jacobi 

Exact 

2 

1 

12.00003 

9.86960 

3 

1 

10.00003 

9.86960 


2 

59.99989 

39.4784 

4 

1 

9.875132 

9.86960 


2 

41.99992 

39.4784 


3 

170.1317 

88.8264 

5 

1 

9.869781 

9.86960 


2 

39.76482 

39.4784 


3 

102.1309 

88.8264 


4 

380.2456 

157.9137 


Table 10 

Nondimensional bending frequencies, X*/^, of a uniform, hinged- 
free beam with N collocation intervals; (X = w ^mRl^l{EI)R) 


N 

Mode 

Jacobi 

Exact 


1 

24.00014 

15.41821 


1 

16.07796 

15.41821 


2 

105.5557 

49.96487 

4 

1 

15.46149 

15.41821 


2 

56.60144 

49.96487 


3 

295.1689 

104.2477 

5 

1 

15.42012 

15.41821 


2 

50.88158 

49.96487 


3 

131.8585 

104.2477 


4 

628.4688 

178.2697 
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Table 11 


Nondimensional bending frequencies, of a uniform, free-free 
beam with N collocation intervals; (X = I{EI)r) 


N 

Mode 

Jacobi 

Exact 

3 

1 

24.49481 

22.37332 

4 

1 

22.61198 

22.37332 


2 

76.68076 

61.68503 

5 

1 

22.38678 

22.37332 


2 

64.21640 

61.68503 


3 

174.3449 

120.9027 


Table 12 


Nondimensional buckling speeds, 0, of a uniform, inward-oriented, 
rotating beam. Jacobi integrating matrix solutions with N colloca- 
tion intervals. (H = a*\/X7) 



WKB 

Direct 

Integration 

Hodges 

Ritz 

FEM 

State Vector Method 
iV = 5 iV = 4 iV = 3 

N = 2 

.2 

.35 

.35 

.35 

.35 

.35 

.65 

.5 

1.13 

1.13 

1.13 

1.13 

1.13 

1.19 

1.0 

2.96 

2.99 

2.99 

2.99 

2.99 

3.05 

1.5 

5.35 

5.38 

5.38 

5.38 

5.38 

5.46 

2.0 

8.27 

8.19 

8.19 

8.19 

8.19 

8.30 

3.0 

14.91 

14.87 

14.87 

14.87 

14.89 

15.05 

4.0 

22.70 

22.77 

22.77 

22.77 

22.81 

23.03 
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Chapter 6 


Divergence and Aeroelastic Lift of 
Composite Wings 


Static aeroelastic behavior, which includes both divergence instability 
and aeroelastic lift distribution, is an important consideration in the design of 
elastically flexible lifting surfaces. An innovative approach to aeroelastic design 
is provided by the concept of aeroelastic tailoring, which addresses the problem 
of designing a flexible lifting surface to take advantage of structural deformation. 
Essentially, one strives to control the deformation, and thus the load distribution, in 
a way that enhances aerodynamic performance. Because of their unique directional 
properties, advanced composite materials prove to be an important ingredient in 
many tailored designs. 

Along with the increased design flexibility allowed by composite materials, there 
comes an additional complexity in the structural analysis. The hybrid state vector 
method discussed in previous chapters provides a simple, yet powerful tool for 
analyzing large aspect ratio composite lifting surfaces. The primary focus of this 
chapter is on applying the anisotropic beam equations developed in Section 2.3 to 
the static aeroelastic analysis of forward swept composite wings. 

The potential benefits of aeroelastic tailoring applied to forward swept com- 
posite wings have been thoroughly examined by Weisshaar [4,5,50]. In order to 
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facilitate the verification of the hybrid state vector solutions, the problems con- 
sidered in this chapter are similar in nature to those discussed by Weisshaar for 
uniform geometry wings. In the following analyses, aerodynamic loads will be cal- 
culated from modified strip theory aerodynamics. 


6.1 Divergence of a Forward Swept Composite Wing 

The cantilevered composite wing structure considered in the following analysis 
is identical to the lifting surface model presented in Fig. 1 of Chapter 2. The 
aeroelastic equations apply to aerodynamic sections taken normal to the swept 
structural reference axis. A detailed derivation of the differential equations for this 
plate-beam model, including important assumptions, is presented in Section 2.3. 
It should be noted that for static aeroelastic stability calculations, these equations 
are, in fact, perturbation equations, with the state vector containing perturbation 
quantities. 

For the sake of simplifying the presentation, it is assumed that transverse shear 
terms can be neglected. Another important assumption to be used in the analysis 
is that the composite wing structure, which is modelled as an equivalent composite 
plate, can be considered as having the same properties as a midplane symmetric 
laminate. Primarily, this means that the 'B*j coupling compliances in Eq. (2.28) will 
be taken as zero. 

For fixed wing problems, it is clear that no appreciable external spanwise loads 
exist. In the event that coupling compliances are included in the analysis, there 
will be induced spanwise loads that arise from satisfying the fifth equation in Eq. 
(2.28). These coupling-induced spanwise loads can be directly calculated if one has 
solutions for the static spanwise midplane deformations «o- If solutions for «o are 
not readily available, the coupling-induced loads can be approximately determined 
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by assuming that the spanwise strains and strain gradients are negligible. 

By following the solution procedure described in Chapter 4, the T matrix 
corresponding to Eq. (2.28) can be easily determined for a cantilevered composite 
wing. When the B*,- terms are assumed to be zero, this matrix can be written as 


T = 


L'A^iLj 


0 

0 

0 


L'Z?tiLt -VDuLf 
-L^'DuLl L^'DuLf 
L'D*M -L'^tsLf 


0 


L'i>33Lt 


( 6 . 1 ) 


A steady aerodynamics matrix based on strip theory can be extracted from the 
unsteady strip theory results in Appendix D by taking s* = 0 , C( 0 ) = 1 . 0 , and 
R = 1 . 0 . The resulting aerodynamic matrix is given by 


rO 0 0 0 1 


Afo 


0 0 0 0 

0 ' 0 'La 

0 ' 0 ' 


(6.2) 


where the terms in Afd are obtained from the corresponding terms of Eqs. (D.1-D.2) 
in Appendix D. For example, X' = '£ 7 ( 5 * == 0 , X), where X is a nondimensional 
dynamic pressure parameter (see Appendix D). For large aspect ratio and moderate 
sweep angles, the terms in Eqs. (D.3-D.4) involving ' L? and ' Mr are negligible. 

Substituting Eqs. ( 6 . 1 - 6 . 2 ) into Eq. (4.29) leads to the divergence eigenvalue 
problem 


■[•Gii 

Gi2‘ 

1 

‘I o' 

in 

. .G21 

G22. 

X 

.0 I. 

Juj 


= 0 , 


(6.3) 
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where 


G„ = -L'UjiLf ' 4 + L':Dl3Lt' 

Gi 2 — — L La + L'Di3Li' Ma 

(6.4) 

G21 = -L'C^sLj^' + L'S^sLi' 

G22 = -L':ot3Lf ' ila + L'P^aLt' . 

Since the divergence problem being considered involves neither axial loading 
nor constitutive terms that couple axial deformations with other degrees of freedom, 
one might expect that the use of simple strip theory aerodynamics, as given above, 
should allow representation of the eigenvalue problem in Eq. (6.3) in terms of an 
effective angle of attack. This is indeed the situation. 

In order to carry out the reduction to effective angle of attack, first note that 
because of the nature of the aerodynamic terms one can write 


Gil 

Gi2 


■— Gi2 

G12' 

' tan A 0 

.G21 

G22. 


— G22 

G22- 

. 0 I. 


Second, note that a constraint equation involving the effective angle of attack Oe 
can be written as 


tte == a — 7 tan A . (6.6) 

This constraint equation expresses the fact that the effective angle of attack for a 
swept wing consists of the actual twist of the wing plus an induced angle of attack 
that arises from lateral bending. From this constraint equation, one can obtain the 
transformation 



Table 13 


Geometric, aerodynamic, and structural parameters for a uniform 
composite wing 


ELASTIC MODULI : 
E 2 JE 1 — 0.1 
G 12 /E 1 = 0.0373 
U 12 — 0.25 


WING PROPERTIES : 
e/bR = 6.67 
bR = 1.0 
1.0 
ao = 27 t 

Oc — 0.5 
a = 0.34 


where 


u-i = 


— ' tanA 
0 


I 

I. 


(6.8) 


By applying this transformation in the form of a similarity transformation to Eq. 
(6.3) while making use of Eq. (6.5), one obtains a reduced eigenvalue problem that 
can be written as 


[G 22 - (' tan A)Gi 2 ]{ae} = (1/X){ae}. (6.9) 

Details on the transformation process used in obtaining Eq. (6.9) are presented in 
Appendix F. Note that either a similarity transformation or a congruence transfor- 
mation must be used since these transformations do not affect the eigenvalues. 

In order to verify solutions given by Eq. (6.9), divergence eigenvalues were 
obtained for uniform, forward swept, composite wings with cantilever boundary 
conditions. Verification was established by comparing the trends of the results 
with analytical solutions given by Weisshaar [4]. The data used in the analysis 
are presented in Table 13, and are taken from a uniform, composite wing analyzed 
by Housner and Stein [51]. These data are also used in a composite wing flutter 
analysis presented in Chapter 7. 
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From the material parameters listed in Table 13, angle-ply laminate stiffnesses 
were calculated for a midplane symmetric composite plate by using the methods 
given in Appendix C. A special form of these calculations for angle-ply laminates is 
given in Section 4.4.4 of Jones [25]. Box beam type structures are easily analyzed 
by carrying out the analyses in terms of an equivalent composite plate. Divergence 
solutions were investigated for laminates having nondimensional stiffnesses equiv- 
alent to an angle-ply laminate with a given number of layers. Since the magnitudes 
of the composite laminate coupling terms vary with the number of layers, this allows 
investigation of the effects of coupling stiffness on divergence velocity. The max- 
imum values of coupling are obtained for a laminate with a single layer (i.e., the 
same as multiple layers with all layers oriented in the same direction). Thus, a single 
layer angle-ply laminate provides a limiting case solution. As the number of layers 
increases, the coupling terms decrease and the laminate tends toward quasi-isotropic 
behavior. 

Fig. 9 displays the variation of the nondimensional bending and torsion com- 
pliances for a single layer angle-ply laminate as the fiber orientation angle varies. 
Material properties were taken from Table 13. The orientation angle 0, positive as 
shown in Fig. C-2, is measured with respect to the laminate reference axis, which 
in the following problems will be assumed the same as the structural reference axis. 
For analyzing a composite wing, the structural reference axis is chosen as the mid- 
chord of the structural box, since the composite stiffness expressions are developed 
with respect to this axis. One should note that and D 33 are symmetric with 
respect to ^ = 0 , whereas the coupling term is antisymmetric in 0. The com- 
pliances in Fig. 9 are nondimensionalized as shown in Eqs. (2.30) and (2.32), with 
the reference stiffness chosen to be Du evaluated at 5 = 0 . 
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FIBER ANGLE, 6 (DEG) 


Fig. 9. Nondimensional compliances for a symmetric angle-ply 
laminate with a single layer 

Presented in Fig. 10 are the nondimensional divergence velocities of a uniform, 
cantilevered, composite wing that is swept forward 30 degrees (A = —30°). Data 
are taken from Table 13. The nondimensional divergence velocities are plotted as a 
function of the fiber orientation angle for the equivalent of a midplane symmetric 
angle-ply laminate. The reference divergence velocity for this figure is the velocity 
for a one layer laminate with 6 = — 9(f . Three different curves are shown cor- 
responding to coupling stiffnesses equivalent to laminates with one, five, and fifteen 
layers. As can be seen, the predicted divergence velocities dramatically increase for a 
range of fiber orientation angles greater than zero. These angles correspond to fibers 
that are oriented ahead of the structural reference axis. Some of the fiber orien- 
tation angles yield negative eigenvalues that give imaginary divergence velocities. 
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FIBER ORIENTATION ANGLE. 6 (DEG) 


Fig. 10. Nondimensional divergence velocities of a uniform, can- 
tilevered, composite wing — symmetric angle-ply with 1, 

5, and 15 layers. (A = -30°, \rbf = X(a=-3O»,0=-9O°) ) 

indicating that divergence does not occur for these fiber angles. The largest range 
of infinite divergence velocities is exhibited by the single layer equivalent laminate, 
which possesses the largest value of bending-twist coupling. Although the bending- 
twist coupling is not the only factor determining the divergence velocity, it is clear 
that it has a dominant efiFect. The physical reason underlying the importance of 
the coupling stiffness is that it can give a lifting surface a washout property, mean- 
ing that bending of the surface, and resulting twist, act to alleviate the excessive 
buildup of aerodynamic loads. 

Similar divergence velocity solutions are given in Fig. 11 for a wing that is 
swept forward 60 degrees (A = —80°). Again, three curves are presented, each 
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FIBER ORIENTATION ANGLE. 6 (DEG) 


Fig. 11. Non dimensional divergence velocities of a uniform, can- 
tilevered, composite wing — symmetric angle-ply with 1, 

5, and 15 layers. (A = — 6Cf , \ref = X(A=-3ff»,0=-flOP) ) 

corresponding to an equivalent angle-ply laminate with a specified number of layers. 
The reference velocity is the same as that in Fig. 10, which is the divergence 
velocity for a single layer laminate having a fiber orientation of ^ = -90" and a 
forward aerodynamic sweep of 30 degrees. For this example, no infinite divergence 
velocities are predicted since the coupling and bending stiffnesses are not large 
enough to completely override the aerodynamic loading. The maximum values 
of divergence velocity still occur at fiber orientation angles slightly ahead of the 
structural reference axis. An interesting point to note from Fig. 11 is that the 
result for a laminate with 15 layers, which has only a small amount of coupling 
stiffness, is nearly symmetric about ^ = 0. This indicates that the behavior of a 
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wing with 60 degrees of forward sweep is also strongly influenced by the bending 
stiffness. The bending stiffness is symmetric in 6 and peaks at ^ = 0. 

For composite materials that are different from the one on which Figs. 10 and 
11 are based, one can expect the locations of the peaks in the divergence velocity 
to shift slightly. These shifts are possible because other composites may exhibit 
different ratios of coupling to bending stiffness. The basic trends demonstrated are 
expected to be valid, however, for any type of unidirectional laminated composite. 
The trends demonstrated in Figs. 10 and 11 are corroborated by the analytical 
results presented by Weisshaar [4]. 


6.2 Aeroelastic Lift of a Forward Swept Composite Wing 

Of great interest in aeroelastic design is the equilibrium lift distribution of a 
flexible wing, especially since aerodynamic efficiency is strongly impacted by the lift 
distribution. The lift distribution of a flexible wing is also closely tied to maneuver 
performance as well as to divergence instabilities. That is, the combined bending 
and twisting of a loaded wing can act either to amplify or to attenuate the loading 
and the bending moments associated with disturbances about a nominal equilibrium 
lift condition. 

Successful design of forward swept wings hinges upon the ability to “tailor” 
the lift distribution. Composite materials offer significant advantages in the design 
and construction of lifting surfaces that are tailored to maintain favorable defor- 
mation patterns, and therefore favorable lift distributions, over a range of flight 
conditions. The following analysis investigates solutions for symmetrical aeroelastic 
lift of forward swept, cantilevered wings. These wings can be described by the same 
composite plate-beam model employed in Section 6.1. 

The approach adopted in solving for the lift distribution involves specifying 


- 112 - 



the pitch attitude (angle of attack) of the untwisted, rigid wing. This provides 
the initial load distribution from which the elastic deformation is determined. The 
total deformation of the wing is the sum of the rigid wing attitude plus the elastic 
deformation component. Correspondingly, the total lift distribution of the wing 
(assuming linear aerodynamics) is the sum of the rigid wing lift plus the elastic lift. 
A similar approach has been used by Diederich and Foss [3] to study the lift of 
metallic swept wings, and by Weisshaar [5] to examine the lift of composite wings. 
Further elaboration on solutions for aeroelastic lift distribution can be found in 
Chapter 8 of Bisplinghoff, Ashley, and Halfman [45]. 

Equations that provide the lift distribution solutions for the aforementioned 
composite plate-beam can be obtained by substituting Eqs. (6. 1-6.2) into Eq. (4.26), 
resulting in 


I-XGii 

-XGi 2 ■ r' 

n TGii 

= X 
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(6.10) 
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where the G terms are given by Eq. (6.4). The attitude of the rigid wing is 
specified by 7r and a^. Similar to the reduction of the dependent variables in the 
divergence problem of Section 6.1, the linear system in Eq. (6.10) can be reduced 
to an effective angle of attack variable, «e, since neither axial loading nor axial 
deformation coupling terms are present. By employing the same relationships and 
the same approach presented in Eqs. (6. 5-6. 8), Eq. (6.10) reduces to 

[I-^{ae} = G{aer) (6.11) 


where 


G = X[G22~( tanA)Gi2j. (6-12) 

To verify these results, lift distribution solutions have been obtained from Eq. 
(6.11) for a uniform, swept-forward composite wing. Material property data for the 
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analysis were taken from Table 13, and aerodynamic loads were calculated with the 
aid of simple strip theory aerodynamics. Results are displayed in Fig. 12 for the 
maximum and minimum values of normalized elastic lift distribution, as 

a function of the nondimensional spanwise coordinate x, A symmetric angle-ply 
composite with a single layer has been assumed (i.e., all fibers are oriented in the 
same direction). Since the lift curve slope is assumed constant, the use of strip 
theory aerodynamics allows the lift distributions to be calculated as the ratio of 
elastic wing deflections (effective angle of attack) to rigid wing deflections. For 
solutions obtained in this manner, it is convenient to specify the rigid wing effective 
angle of attack as a unit value. Also specified prior to the solution of Eq. (6.11) is 
the nondimensional dynamic pressure parameter X. The value for this parameter 
was arbitrarily picked to be one half of the divergence dynamic pressure parameter 
for a wing having 60 degrees of forward sweep and an angle-ply fiber orientation 
angle of ^ = — 9(f. 

For the wing used in the analysis, it is assumed that all fibers are oriented at 
an angle 0 with respect to the structural reference axis. The results display for each 
constant wing sweep angle the fiber orientation angle that corresponds to either the 
maximum load amplification or attenuation that can be achieved by orienting the 
reinforcing fibers. Hence, these solutions are an indicator of the maximum amount 
of “tailoring” that can be achieved for a given configuration and material. The 
designs that fall below the dashed rigid wing reference line in the figure are referred 
to as load attenuating designs, while those above the reference line correspond to 
load amplifying designs. It is interesting to note that the maximum amount of lift 
attenuation, which also corresponds to maximum divergence speed, occurs when 
the fibers are oriented ahead of the structural reference axis. Those swept-forward 
wings whose lift distributions fall below the dashed rigid wing reference line are 
in fact displaying the load alleviating property of an isotropic swept-back wing. 




^ = X It 


Fig. 12. Limiting elastic lift distributions for a uniform com- 
posite wing — symmetric angle-ply with a single layer 

These results are verified by an analytical solution developed by Weisshaar [5] for 
normalized elastic lift distribution of uniform composite wings. A comparison of 
the analytical solutions with the approximate Jacobi integrating matrix solutions 
(based on five discretization intervals) shows a maximum deviation from the true 
lift distribution curves of less than one-tenth of one percent. 

Another conclusion that comes from Fig. 12 is that the wing of this example, 
with 60 degrees of forward sweep, never completely reaches the load attenuating 
region. This corresponds to the fact that for A = — 6(f , the wing will always possess 
a finite divergence velocity. For the other wing sweep angles given in the figure, 
there exists a fiber orientation angle for which the result is coincident with the rigid 
wing reference solution. Such designs are referred to as aero-isoelinic and have 

- 115 - 


I 



neither a load attenuating nor load amplifying character. 
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Chapter 7 


Flutter of Isotropic and Composite Wings 


Dynamic aeroelastic behavior will now be examined for both isotropic 
and composite wings having cantilever boundary conditions. The solutions to be 
pursued in the following problems will yield the airspeeds associated with flutter 
instabilities. For composite lifting surfaces, it will be shown that the fiber orientation 
angles have a very strong influence on the points of instability, as well as on 
the subcritical dynamic response. It should be pointed out that the following 
investigations are meant only to demonstrate the solution capabilities available with 
the hybrid state vector approach. A detailed investigation of various types of flutter 
phenomena is not attempted here since the scope of such a study falls outside the 
objectives of the present work. 

In the sample flutter calculations to be presented, a single formulation of the 
flutter equations will be used that is appropriate for anisotropic structures; isotropic 
structures are simply considered as a subset of the anisotropic case. The approach 
used to obtain the subcritical dynamic response and the flutter points involves 
tracing out the complex roots loci of the matrix flutter equations. A description of 
this solution process will be presented in the course of the following analyses. 
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7.1 Flutter Equations 


The general form of the flutter eigenvalue equation to be used was presented 
in Eq. (4.30). It appears as 


I + 


+ ^)] 



(7.1) 


where is the discrete mass matrix for the lifting surface, Cfd is a damping 
matrix, and ^fd is an unsteady aerodynamic matrix. Note that Eq. (7.1) is actually 
the Laplace transform with respect to time of the homogeneous, unsteady aeroelastic 
equations. It is assumed that the equations have been nondimensionalized such that 
s* is a nondimensional, complex-valued, Laplace transform variable given by 



and y£, is a nondimensional complex eigenvector. Also note that in the present 
application X is a dimensionless dynamic pressure parameter whose definition is 
given in Appendix D. 

Precise forms of the matrix terms in Eq. (7.1) are taken to correspond to the 
anisotropic plate-beam equations presented in Eq. (2.28). The unsteady aerody- 
namic matrix, which will be calculated from modified strip theory, is given in 
Appendix D. The T matrix has already been presented in Eq. (6.1) for cantilever 
boundary conditions. It will also be assumed for convenience that the damping 
terms are zero for the following analyses. In nondimensional form, the mass matrix 
in Eq. (7.1) can be expressed as 


M 


FD 


ro 0 0 0 

0 0 0 0 

0 0 Uliunj “^wot 


0 0 ' m,. 


^OtOC'J 


(7.3) 
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where 



( 7 . 4 ) 


The term Xa represents a section mass static unbalance about the structural refer- 
ence axis, and Ta represents a section dimensionless radius of gyration. 

Taking into consideration the zero columns in the mass and aerodynamic 
matrices, it is readily shown for the case of no damping that the flutter eigenvalue 
problem can be reduced to 


(I + G22) 

G23 

G24 

n 

G32 

(I + G33) 
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. G42 
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(I -1- G44)_ 

la) 


where 


G = T[Mfds*^ - <^fd{s*, X )] . ( 7 . 6 ) 

Even further, the flrst row and first column of Eq. (7.5) can be dropped if the lifting 
surface has no aerodynamic sweep. Eq. (7.5) can be numerically evaluated for the 
roots loci of the aeroelastic modes as a function of the dynamic pressure parameter 
X, with instability being associated with complex roots s* having positive real parts. 

The procedure used in the present analyses to trace out the complex roots loci 
was based on a determinant iteration scheme. A determinant iteration method has 
the advantage of being simple to implement, and furthermore, it is not restricted 
to solving a particular form of the flutter equations. If desired, the determinant 


- 119 - 



method can easily solve the flutter equations when the unsteady aerodynamic terms 
are calculated directly from the Bessel function representation of the generalized 
Theodorsen function. 

Since the roots of the flutter equation correspond to the zeroes of the deter- 
minant of the complex- valued matrix in Eq. (7.5), the determinant method simply 
searches for these zeroes. For a given dynamic pressure and a trial root, the deter- 
minant is numerically evaluated from the product of the diagonal terms of the tri- 
angular factor that is obtained by applying Gaussian elimination. Muller’s method 
is used to iteratively search for the complex roots of the determinant. 

A computer code for determinant iteration was developed directly from routines 
available in Chapters 2 and 3 of Conte and de Boor [38]. The resulting program 
worked remarkably well for tracing out the roots loci. Solutions could be started 
by specifying trial roots along the positive imaginary axis and letting the solutions 
converge to the free vibration eigenvalues. By incrementing the dynamic pressure 
and using the roots at the previous dynamic pressure point as iteration starting 
values, a specified number of branches of the complex roots loci could be simul- 
taneously traced out. Any complex valued roots had their conjugates added to the 
function deflation in the Muller routine so that only the upper half-plane of the 
symmetric roots loci was extracted. Zero frequency static divergence roots could 
also be extracted with determinant iteration. 

One point to be aware of is that small sized flutter equations (i.e., not too many 
degrees of freedom) will usually not encounter problems with determinant evaluation 
since the range of determinant values tends to be reasonable. For larger problems, 
one may find it necessary to employ scaled arithmetic in the determinant evaluation 
in order to prevent overflow or underflow during the computation. Another alterna- 
tive is to keep the size of the flutter equations small by applying modal superposition 
to Eq. (7.5). As a result of standard superposition procedures, the flutter matrix 
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Table 14 


Flutter velocities for an unswept, uniform, 
isotropic wing 




Flutter 

Solution 

Description 

Velocity 

km/hr 

Ref. [52,53] 

Exact Analysis 

494 

Ref. [51] 

25 finite— difference points 

483 

State Vector 

Jacobi I.M., 2 intervals 

486 

State Vector 

Jacobi I.M., 3 intervals 

494 

State Vector 

Jacobi I.M., 4 intervals 

494 


will be formed by premultiplying the matrix in Eq. (7.5) by the transpose of the 
modal matrix and postmultiplying it by the untransposed modal matrix. 


7.2 Flutter of an Isotropic Wing 

Flutter calculations using the state vector approach were verified through 
comparisons with analytical solutions presented by Goland [52] for an unswept, 
cantilevered, isotropic wing. (Corrections for Ref. [52] are given in Ref. [53].) The 
appropriate data for this problem were taken from Ref. 44. Table 14 gives a 
comparison of numerical flutter results from Eq. (7.5) with those of Goland and with 
a numerical solution of Housner and Stein [5lj. The state vector solutions are based 
on Jacobi integrating matrices corresponding to two, three, and four collocation 
intervals, and unsteady aerodynamic calculations employ a rational approximation 
of the Theodorsen function due to R.T. Jones. Details on the unsteady aerodynamic 
strip theory can be found in Appendix D. 

It is readily seen that the hybrid state vector solutions are quite accurate 
and converge rapidly to the exact solution as the number of collocation intervals 
increases. Similar solutions based on Newton integrating matrices also yield highly 
accurate results that are virtually identical with the Jacobi solutions. The flutter 
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Table 15 


Inertia parameters for the uniform composite wing of Table 13 

Xa ~ — 1*0 

ra — 0.29 f^wa — 0.03 

H = 8.0 maa — 0.0065 


solutions for this particular problem depend heavily on the accuracies of the lower 
frequency coupled bending-torsion modes of free vibration. The high accuracy with 
which the flutter velocities are determined demonstrates the fact that the integrating 
matrix solutions contain precise information for these lower modes. 

7.3 Flutter of a Composite Wing 

Flutter solutions will now be obtained for an unswept, uniform composite wing. 
The wing geometric, aerodynamic, and structural data to be used in the flutter 
solutions can be found in Table 13 and Fig. 9 of Section 6.1. (Fig. 9 contains data 
for a single layer laminate only.) Inertia data are given in Table 15. This data 
corresponds to a wing that was analyzed by Housner and Stein [51]. The analysis 
to be given here, however, includes the effect of bending-twist coupling represented 
by the Djj term. It is assumed that the composite layup is equivalent to a midplane 
symmetric angle-ply laminate, for which the effective coupling decreases with an 
increasing number of layers. As in the preceding section, the unsteady aerodynamic 
loads are calculated from the modified strip theory presented in Appendix D, and 
the Theodorsen function is calculated from the R.T. Jones rational approximation. 
No compressibility corrections are applied to the unsteady aerodynamic terms. 

Plotted in Fig. 13 is the nondimensional flutter dynamic pressure parameter 
for the unswept, uniform composite wing as a function of the fiber orientation 
angle 0. The results are normalized by the flutter dynamic pressure at ^ = —90°. 
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FIBER ORIENTATION ANGLE, d (DEG) 


Fig. 13. Effect of fiber orientation on flutter dynamic pressure 
for an unswept, uniform composite wing — symmetric 
angle-ply laminate with 1, 3, 5, and 15 layers. {\ref = 

^(0=-9O“)) 

Several curves are displayed corresponding to symmetric angle-ply laminates having 
1, 3, 5, and 15 layers. The single layer equivalent laminate possesses the maximum 
amount of bending twist coupling, while the 15 layer equivalent possesses very little 
coupling and exhibits essentially quasi-isotropic behavior. It is apparent from Fig. 
13 that in the limit of quasi-isotropic behavior, the flutter speed for the unswept 
wing follows the torsional stiffness variation, with a peak in the flutter dynamic 
pressure occurring at maximum torsional stiffness corresponding to (? = 45°. These 
nondimensional results for the 15 layer laminate compare very closely with the 
quasi- isotropic results given in Fig. 9 of Housner and Stein [51]. 
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In contrast with the quasi- isotropic limit, the behavior of a single layer equiv- 
alent symmetric angle-ply laminate demonstrates a limiting case solution involving 
maximum bending-twist coupling. The asymmetry of the bending-twist coupling 
term about ^ = (f is quite apparent from Fig. 13. The sharp peak and the rapid 
changes in the flutter dynamic pressure for this case are caused by the important 
participation of the coupling term in the solution. The rapid variation and the 
sharpness of the peak are quickly diminished when more layers are taken. Since 
practical composite layups tend more toward quasi-isotropic behavior, it is not an- 
ticipated that this single layer limiting case solution will necessarily be attained in 
practice. 

The flutter curves presented in Fig. 13 were obtained directly from Eq. (7.5) by 
using the previously described root tracking scheme based on determinant iteration. 
Both Jacobi and Newton integrating matrices for two and three collocation intervals 
were employed in the flutter calculations. By using the two different interval sizes, 
the convergence of the flutter solutions as a function of the discretization level 
could be checked. For the composite wing flutter problem examined here, it was 
concluded that three intervals gave sufficient convergence, as it did for the case of 
isotropic wing flutter in Section 7.2. Because of the small number of discretization 
intervals needed, both the Jacobi and Newton integrating matrix solutions gave the 
same result. It should also be noted that since solutions were obtained directly from 
Eq. (7.5), modal superposition was not employed and the question of the number of 
modes used in the analysis does not arise. For flutter problems in which much more 
discretization is needed to adequately describe the structure, modal techniques can 
be applied to keep the flutter solutions within manageable proportions. 

To obtain each flutter dynamic pressure, the complex roots loci (s* = cr*+ju!*) 
for that value of fiber orientation angle were traced out beginning with zero dynamic 
pressure and continuing until one of the branches had roots with positive real parts. 
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thus indicating a dynamic instability. Figs. 14 and 15 give typical examples of 
roots loci for two different sets of wing parameters. Fig. 14 gives the roots loci 
corresponding to the highest peak in the flutter dynamic pressure for the three 
layer laminate solution plotted in Fig. 13. This peak occurs at a fiber angle of 
0 = l(f. Fig. 15 gives the roots loci for a single layer equivalent laminate when 
the mass ratio is taken to he n — 64 and the fibers are oriented at ^ == 45°. For 
both of these figures only the three lowest root branches are shown. These two 
figures are indicative of two different types of instability. The instability in Fig. 14 
is approached rather slowly as dynamic pressure increases, whereas the instability 
shown in Fig. 15 is approached very rapidly with a small change in dynamic pressure 
and involves strong coupling between two of the wing modes. 

By tracing each root branch starting with zero dynamic pressure, it was possible 
to avoid the program logic necessary to detect discontinuities in the flutter speed 
as 6 varies. Such discontinuities were encountered by Housner and Stein [51] for 
certain values of the mass ratio parameter n, but not for /i = 8, which was used in 
obtaining Fig. 13. With the implementation of such detection logic, flutter solutions 
as a function of 6 can be carried out by using the flutter dynamic pressure at an 
adjacent value of as a starting point for the next flutter solution. 

The current investigation did not attempt to make a complete examination 
of the instabilities associated with composite wings. A more thorough study of 
the instability boundary must involve simultaneous consideration of divergence 
and flutter and should include the effects of rigid body modes. As indicated by 
Weisshaar [50], the stability boundary of cantilevered, forward swept composite 
wings will be determined for a wide range of fiber angles by the low divergence 
velocities associated with those fiber orientations. Also indicated in Ref. [50] are 
various changes in the mode of flutter instability as the fiber orientation changes. 
The presence of these different modes of instability as the wing parameters vary 
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is worthy of further investigation. It is felt that the solution methods presented 
here will provide a convenient tool for conducting further research on instability 
phenomena of composite lifting surfaces. 
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15 



Fig. 15. Roots loci of aeroelastic modes for a uniform composite 
wing. (Single layer laminate; fi = Qi, 0 = 45°) 
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Chapter 8 


Concluding Remarks and Recommendations 


The previous chapters describe a hybrid state vector method for solving 
the differential equations arising in structural dynamics and aeroelasticity. The 
method proves to be very versatile and can be applied to the solution of any form 
of ordinary differential equation. 

Solutions can be obtained in a consistent fashion for any size problem by 
working with the state vector form of the differential equations. An integrating 
matrix method is used to discretize the differential equations, yielding standard 
linear matrix equations from which one obtains the desired solutions. It is also 
shown that for simple problems, the integrating matrix can be applied directly to 
other forms of the differential equations. 

It is demonstrated that a convenient form of the state vector equations can be 
obtained from a variational formulation of the structural equilibrium equations. The 
equations given by this formulation have properties that are useful for numerical 
calculations. For structural problems, the state vector equations can be partitioned 
corresponding to generalized force and displacement variables. By applying matrix 
partitioning techniques, the solution state vector can be reduced to the displacement 
state variables only. This reduction process can often be carried out analytically 
to yield expressions for the direct calculation of the terms in the reduced matrix 
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equations. For problems in which it might not be convenient to carry out analytical 
reductions of the equations, the calculations can be performed numerically. 

As demonstrated in Chapters 5, 6, and 7, the hybrid state vector method is 
applicable to a variety of problems and boundary conditions. Results for these 
problems indicate that high accuracy is easily achievable. By employing high 
order polynomial approximations, the integrating matrix solutions are capable of 
providing sufficient accuracy with minimal discretization. The matrix operations 
required by this solution procedure are straight forward, easily programmable, and 
allow for efficient problem solution. 

A theory of integrating matrices is presented and a calculation procedure is de- 
veloped for maximum precision integrating matrices that are based on orthogonal 
polynomials. Discussions are given for several types of integrating matrices and 
some of their properties. Much flexibility is available in the types of integrating 
matrices that can be derived and used for special applications. Each type of in- 
tegrating matrix has its own uniques properties that make it more or less applicable 
to a particular situation. For the problems examined in this work, the integrat- 
ing matrices based on Jacobi polynomials have proven to have good convergence 
properties and are capable of very high accuracy solutions. Newton based integrat- 
ing matrices, which are convenient for applications requiring evenly spaced grid 
points, also provide very good accuracies. The convergence of solutions based on 
Newton matrices, however, tends to be oscillatory in character unless the number 
of grid points is somewhat larger than the minimum number of points required by 
a given order of approximating polynomial. Tabulations are presented in Appendix 
B for both Newton and Jacobi integrating matrices. 

Additional work needs to be devoted toward extending the theory of hybrid 
state vector solutions to two and three dimensional problems. It should prove 
to be convenient in the development of multidimensional integrating matrices to 
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introduce the familiar concept of multidimensional shape functions. This concept 
is routinely employed in finite element analyses. For irregular shaped regions, it is 
anticipated that a transformation to a simpler domain will be useful in performing 
the integrations. That is, integrating matrices are more easily developed for these 
simpler domains. Again, this concept is very much similar to practices currently 
in use for isoparametric finite element solutions. It is also possible to extend the 
present method to two dimensional problems by applying separation of variables, 
thus allowing reduction of the partial differential equations to ordinary differential 
equations. For multidimensional problems in general, it should prove worthwhile to 
investigate the use of integration methods that are especially suited to multidimen- 
sional integration. 

Other types of integrating matrices should be investigated for special types of 
problems. For instance, semi-infinite boundary value problems should easily yield to 
numerical solution by integrating matrix. In this connection, it would be worthwhile 
to examine the prospect of developing integrating matrices from the orthogonal 
polynomials that are normally used for quadratures on semi-infinite domains. 

Further studies are also warranted for applications of the hybrid state vector 
method to nonlinear problems. For materially nonlinear structures in particular, it 
appears that the Reissner variational formulation (written in terms of the hybrid 
state vector) may offer important advantages for numerical solution by integrating 
matrices and is deserving of a more thorough investigation. In the presence of 
nonlinearities, the hybrid state vector formulation, in conjunction with integrating 
matrices, provides a very convenient method for obtaining a compact set of non- 
linear algebraic equations that describe the solution of the nonlinear problem. 

Finally, there are many questions yet to be answered concerning the fiutter of 
composite wings. Some of the more important questions are mentioned in Chapter 
7. The present examination does not attempt to address all of these problems, but 
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rather forms the foundation for further parameter studies of wing flutter. Most 
certainly, the answers to some of these problems will be very useful in preliminary 
aeroelastic design. 
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Appendix A 


Weighting Matrices and Shape Functions 
for Jacobi Polynomials 


The following presents a summary of the calculation procedures for Jacobi 
integrating matrices. Jacobi integrating matrices are developed from Lagrange 
interpolating polynomials that can be written in terms of the normalized Jacobi 
polynomials pn Included in this formulation is a detailed description of the 
integration methods required to calculate the Jacobi weighting matrices. In addi- 
tion, the interpolating polynomials are presented in the form of shape functions 
(basis functions). Discussions on shape functions can be found in Zienkiewicz [54] 
and Gallagher [55]. 

First, a distinction must be made between the notation in Chapter 3 and 
the notation appearing in this appendix. In Chapter 3, n denotes the number of 
subintervals in the interval of integration, with n + 1 being the total number of grid 
points on the interval. In this appendix, however, n refers to the number of grid 
points in the interior of the interval of integration (i.e., excluding end points). This 
notational change prevents undesirable complication of subscripts and limits and 
allows the use of an accepted notation for Jacobi polynomials. Since the notational 
change is confined entirely to the calculations in this appendix, no confusion should 
arise. It should be noted that whenever n is used as a subscript in a polynomial, 
it denotes the degree of the polynomial; however, the degree is the same as the 
number of interior grid points. 

All of the integrations and interpolations will be written for the normalized 
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interval [—1,1], and the subscript k will designate the quadrature points on this 
interval. The quadrature points Xh are fixed grid point locations that are determined 
by the zeroes of the appropriate Jacobi polynomial. These quadrature points are the 
same as the abscissas for Lobatto integration, which can be found, for instance, on 
page 920 of Abramowitz and Stegun [40]. Appendix B lists these abscissas for each 
Jacobi integrating matrix. Note too, that a distinction will be made between the 
unnormalized Jacobi polynomials and the normalized polynomials 

These polynomials are a special case of the general Jacobi polynomial 

= \/^ (A.l) 


where the normalizing factor 8n is 


^ (a + + 2n + l)n!r(a + /? + »» + 1) 

” 2“+/?+ir(a + 71 + i)r(yi? + 71 + i)~ ■ 

For a = ^ = I, the normalizing factor can be reduced to 
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_ (2n + 3)(n + 2) 
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(A3) 


Some more preliminaries on Jacobi polynomials are necessary before discussing 
shape functions and integrating matrices. First, the Jacobi polynomials 
can be calculated recursively via the formula 


Pn = (n + l)ro(a:) 


(A.4) 


with 


'•m-lCa;) = 1 + 


(n + 1 — m)(n + 2 + m) 
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The recursive calculation is repeated for m = n, n — 1, . . . , 2, 1, beginning with 


r„{x) == 1. By means of Eqs. (A.4-A.5) one can show that 



(®) = 1 


(A.6) 


and 

= (n + 1). (A.7) 

Using identities listed on page 777 of Abramowitz and Stegun [40], one also finds 
that 

= (-!)"(« + 1), (A8) 

and from the differential relations given on page 783 of Ref. [40], a useful derivative 
term is expressible in the form 




n + 1 




(A.9) 


From the recursion relations given in Eqs. (A.4-A.5), the following series expansion 
of was developed for use in integrations: 


= (n + 1)1 


1+ x: *(*-«' 

m=l 


(AlO) 


where 
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(A.ll) 


Note that the summation on the right-hand side of Eq. (A. 10) vanishes whenever 
n < 1. One final definition that’s needed in the following discussion is the definition 
for the leading coefficient of an orthonormal polynomial. For the normalized Jacobi 
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polynomials with a = ^ = \, the leading coefficient, an (i.e., the coefficient of the 
highest order term a:” ) can be expressed as 


On 


rr ( 2 » + 2 )! 

^ ” 2”n!(n + 2)!' 


(A12) 


Turning now to interpolation, it is assumed that a sufficiently smooth function 
f{x) can be reasonably approximated by a Lagrange interpolating polynomial. If 
the end points of the interval [—1,1] are included in the interpolation, then an 
approximation to f{x) can be expressed as 


fix) ^ 


- 1 ) 


fixk) + 


Pn’^\x)[x^ - 1 ) 


k=\ [x - Xk)p\^’^^'{xk)ixl - 1) y=l (a; - dj)p^rt’^\dj){2dj) 


-/(dy) (A13) 


where the end points of the interval are d\ = — 1 and d 2 = 1. This particular 
form of the Lagrange interpolating polynomial is convenient for use with orthogonal 
functions. 

The interpolation given by Eq. (A. 13) can be equivalently written in terms of 
shape functions, JJ{x), as 


fix) ^ Mkix) M+iix)\{f) (A14) 

where 

(/} = (/(-l) /(+!))’■. (A15) 

If interpolation is to be performed for multiple points, then the left hand side of 
Eq. (A. 14) becomes a column vector and the shape function terms are written in 
the form of a rectangular matrix rather than a row matrix. 

The shape functions, which provide a convenient way of expressing Eq. (A. 13), 
can be simplified by applying Eqs. (A.6-A.8). As a result of this simplification, the 
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end point shape functions appear as 


and 






For internal points Xf. on the interval [—1, 1] the shape functions are 


= 




(A16) 


(A17) 


(A18) 


where the weighting function i?(a:) = a;^ — 1. This weighting function is the same 
as mentioned in Section 3.1. Note too, that the prime in the denominator of Eq. 
(A. 18) indicates differentiation with respect to x. To facilitate further manipulation 
of the shape functions, the superscript (1, 1) will be dropped in the remainder of this 
appendix. That is, p« and p|/’*^(ar) will now be referred to simply as pn{x) and 
Pn{x). 

Unfortunately, Eq. (A. 18) is very inconvenient for numerical calculations. The 
Christoffel-Darboux identity, however, proves useful in reducing the interior point 
shape functions to a much more desirable form (see, for example, Krylov [39], p. 
103). By using the fact that Pn{xk) = 0 when evaluated at its zero points, x^, a 
useful form of the Christoffel-Darboux identity can be written as 

x: pA-)p.M - 

g_0 <*n+l ® X/i 

If Eq. (A. 19) is first multiplied through by t?(ar) = a:^ — 1, then some rearranging 
gives 

pn{x)d{x) _ anBk{x) 
x-xk a„_iF„_i(x^) ’ 
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where 


n— 1 

h{x) = Y1 ^8/"e(afjfcMa:)F*(ar) . (A21) 

e=0 

To write Eq. (A.20) in its present form, use was made of the standard recurrence 
relation for orthonormal polynomials, which allows one to obtain the equality 


~Qn+l 

O'nPn+xi.^k) ®n— 1-E l(*A;) 


(A.22) 


By substituting Eq. (A.20) into Eq. (A. 18), one finally obtains a convenient 
form for the shape functions 




-8^fe(a:) 

(n + l)(n + 2){Pn-\{xk)f 


(A23) 


where Eq. (A.3), Eq. (A.9), and Eq. (A.12) have been employed to arrive at Eq. 
(A.23). The final form of the shape functions are represented by Eq. (A.16), Eq. 
(A. 17), and Eq. (A.23). With the aid of the recursion relations given by Eqs. (A.4- 
A.5), the shape functions, used in conjunction with Eq. (A. 14), provide a convenient 
way of numerically performing interpolations. These interpolation expressions also 
provide the foundation for calculation of the Jacobi weighting matrices. 

It is useful to note that the interpolations can be carried out in matrix notation 
after defining a matrix of shape function values. For example, if one chooses a 
number of fixed points x = xj, {j — 1,,. ,,t) at which to interpolate, then the tXn 
matrix for the interior point shape functions, is given by 

‘ (^-24) 

where the t X n matrix is written as 

\^j,^\^«Pe,k\' (A.25) 
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The matrix product on the right-hand side of Eq. (A.25) performs the same sum- 
mation as the right-hand side of Eq. (A.21). is a f X n matrix given by 




(A.26) 


and is an n X n matrix calculated as 

= l«*F«(a:;t)]. (A.27) 

In addition to f-Vy,*], the following column vectors arise from evaluating the end 
point shape functions via Eqs. (A.16-A.17): 


= {>/_l(a:y)}, (A28) 

{>/+y} = {-V+i(xy)}. (A.29) 


Finally, by collecting Eqs. (A.24), (A.28), and (A.29), one obtains the complete shape 
function matrix in the form 


(A.30) 

This matrix of shape function values can now be used for function interpolation in 
the manner indicated by Eq. (A. 14). 

Having obtained expressions for the interpolation of f{x), the next step is to 
calculate the weighting matrices from these approximations. But first, a convenient 
notation must be agreed upon for the subintervals over which the integrations will 
be performed. In Section 3.1, it was convenient to label subintervals by grid points. 
That is, a particular subinterval was referred to by [ar^, ar,+i], where the subscripts t 
and 14-1 designated the two consecutive grid points on that subinterval. In discussing 
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the integrations required for Jacobi weighting matrices, however, it will greatly 
simplify the notation to refer to the same subinterval by [ly, ary^] (y = 1, , . . , n + 1), 
where j denotes a particular interval rather than a grid point. 

The elements in a Jacobi weighting matrix can be determined by integrating 
the approximate expression for f{x) given in Eq. (A. 14). These integrations, which 
will be carried out over subintervals, are written in the general form 


I f{x)dx= I g{x)'d{x) dx . (A.31) 

Jxj Jxj 

As mentioned in Chapter 3, a maximum precision quadrature will result only if g{x) 
is orthogonal to t?(a;) over the total interval of integration. The total interval of 
integration here is the normalized interval (—1, 1]. Since the weighting function that 
arises during interpolation is ^(x) — x^ — 1, the Jacobi polynomials are the proper 
choice for the orthogonal function. 

Equation (A. 14) is now substituted into Eq. (A.31). When the shape function 
definitions in Eq. (A.16), (A.17), and (A.23) are taken into account, one can define 
the integrals 


and 




J^-i(x) dx — 


(-ir‘ 

2(n + l)/x, ' 


l)Pn{x)dx 


(A.32) 


and finally, 

Cjh — I Ni-{x)dx = 

J X 2 


-8 


(n + l)(n + 2)(P„_i(a!,t)) ^ 


/ Bk{x)dx. 
J X: 


(A.33) 


(A34) 


In order to isolate the integral portion of the terms in Eqs. (A.32-A.33), it is 
convenient to rewrite these equations as 


(-!)-» 

2(« + l)^-^ 


(A.35) 
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and 


2(n + l)^+^’ 


(A.36) 


where 


and 



(A.37) 


f j+ 

Q+j= I {x + l)Pn{x)dx. (A.38) 

JXj- 

With the aid of Eq. (A.21), it is possible to recast Eq. (A.34) in a similar manner to 
obtain 


where 


and 




(n + l)(n + 2)(P„_i(a:^)f 


n-1 

= E 

e=0 


(A.39) 


(A.40) 



(A41) 


The next step in calculating the weighting matrices requires the development 
of expressions for the integrals in Eqs. (A.37), (A.38), and (A.41). Appropriate 
expressions for these integrals are derived with the help of Eqs. (A.lO-A.ll), which 
will allow term by term integration. Therefore, substitution of Eq. (A. 10) for Fn(a;) 
in Eq. (A.37) yields the integration 


Q-j = (n + 



— l)dx + 



{x-l)^+^dx . 


(A.42) 
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Carrying out this integration provides the result 


2-y - (»+i) i(x|+ - ^ c{^i(«,-+ - ir+" - (»y - ir+“i} 

(A43) 

A similar operation applied to Eq. (A.38) yields the integrals 


f {x+l)dx+ ^ eJJ, f [(* — 1)"* + a:(x — 1)”*] da: 
.3 m=l 


Q+J = (n + 1) 


with the final result of the integration being 


(A44) 


Q+j = (w + 1) 


m=l ^ 


(A.45) 

Proceeding with a similar approach for Eq. (A.41) gives 


^y,8 ==(« + !) 


j2 c +*(>’- ir+'i 

.3 m=l '^^3 


dx 


(A46) 


for which the integrated result is 


m=l 


(A47) 

It is important to note that 0<s<ti — 1 as required by Eq. (A.21); however, the 
summation appearing on the right-hand side of Eq. (A.47) vanishes whenever a < 1. 

The Jacobi weighting matrices can now be constructed from the values for D-j, 
D+j, and Cj^k- Numerical calculation of these values makes use of Eqs. (A.35-A.36) 
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and Eqs. (A.39-A.40) in conjunction with Eqs. (A. 43), (A. 45), and (A. 47). Note that 
and D+j actually form column vectors with a dimension equal to the number of 
discretization subintervals. Similarly, the elements Cj^k form a rectangular matrix 
with a row dimension equal to the number of discretization subintervals and a 
column dimension equal to the number of internal grid points, where the number of 
internal grid points is one less than the number of subintervals. Analogous to Eqs. 
(A.24-A.25), one finds that the matrix can be computed from 

^n-l(*fc)] (A.48) 

where the (n + 1) X n matrix [By j;.] 

\^lk\ = \^U[6sPs,k]- (A.49) 


The second matrix on the right of Eq. (A.49) is the same as Eq. (A.27). And finally, 
with these definitions, the Jacobi weighting matrix is constructed in the form 


'M'n+l 


■ 0 ... 0 ■ 
.{Z)_y} [Cy,,l {B+y}. 


(A.50) 


where, as mentioned in Section 3.1, the first row contains only zeroes. Considering 
the definition of n in this appendix, the dimension of IVn+i is (n + 2) X (n + 2). 
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Appendix B 


Tables of Integrating Matrices 


The integrating matrix L was defined by Eq. (3.11) as 

L = SWn. 

This definition applies to all of the integrating matrices given in this appendix. 

Jacobi Integrating Matrices 

The following Jacobi weighting matrices are for the normalized interval [—1, 1], 
and can be transformed to the interval [0,1] by multiplying by the factor 1/2. 
The coordinate values x, which denote grid points on [—1,1], are listed above the 
horizontal bar in each matrix. Equivalent grid points t on [0,1] are given by the 
linear transformation t = 0.5(1 + a:). For the Jacobi weighting matrices, the subscript 
n refers to the number of discretization intervals, with the integration being exact 
for all polynomials of degree <2n — 1. 


“M/2 = 


•-1. 000 00000 .00000000 l.OOOOOOOQ- 

.00000000 .00000000 .00000000 

.41666667 .66666666 -.08333333 

. -.08333333 .66666666 .41666667. 


{B.l) 


Ws 


•- 1.00000000 -.4472 1360 .44721360 l.OOOOOOOQ - 

.00000000 .00000000 .00000000 .00000000 

.22060113 .37939888 -.06781473 .02060113 

-.07453560 .52174919 .52174919 -.07453560 

. .02060113 -.06781473 .37939886 .22060113- 


{B.2) 
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r - 1.00000000 

-.65465367 

.00000000 

.65465367 

i.oooqoopo- 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.13545686 

.23948954 

-.04347144 

.02127165 

-.00740028 

-.05420686 

.36687883 

.39902700 

-.08319557 

.02615028 

.02615028 

-.08319557 

.39902700 

.36687883 

-.05420686 

. -.00740028 

.02127165 

-.04347144 

.23948954 

. 13545686 - 


r - 1.00000000 

-.76505532 

-.28523152 

.28523152 

.76505532 

1 . 00000000-1 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.09135961 

.16373565 

-.02974921 

.01525535 

-.00894356 

.00328685 

-.03954284 

.26394600 

.29767068 

-.06326350 

.03255895 

-.01154547 

.02310851 

-.07282206 

.33494505 

.33494505 

-.07282206 

.02310851 

-.01154547 

.03255895 

-.06326350 

.29767068 

.26394600 

-.03954284 

- .00328685 

-.00894356 

.01525535 

-.02974921 

.16373565 

. 09135961 - 


( B . 4 ) 


r - 1.00000000 




.00000000 

.06569253 

-.02968808 

.01905508 

-.01164044 

.00587650 

-.00167654 


-.83022390 -.46884879 
.00000000 .00000000 
.11864577 -.02153719 
.19675645 .22624682 

-.05982575 .27026169 

.03261326 -.06238108 
-.01579792 .02613302 

.00443419 -.00697786 


■0000 0000 .46884879 
.00000000 .00000000 
.01119518 -.00697786 
.04815170 .02613302 
.28076604 -.06238108 
.28076604 .27026169 
.04815170 .22624682 
.01119518 -.02153719 


. 83022390 _ 1 . 00000000 - 
.00000000 .00000000 
.00443419 -.00167654 
-.01579792 .00587650 

.03261326 -.01164044 
-.05982575 .01905508 

.19675645 -.02968808 
.11864577 . 06569253 . 


( B . 5 ) 
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- 1.00000000 

.00000000 

.04047503 

-.02205750 

.01555127 

—.01050027 

.00051400 

—.00330210 

.00004205 


-.87174016 

.00000000 

.08078532 

.15151413 

-.04873011 

.02036150 

-.01742762 

.00886234 

-.00245244 


.60170018 —.20020022 


.20020022 


.50170018 


.87174015 1 . 00000000 ' 


.00000000 

.01628154 

.17827881 

.21808803 

-.05570463 

.02840036 

.01335611 

.00360077 


.00000000 

.00850816 

.03746083 

.23145357 

.24624152 

-.05153854 

.02067231 

.00540841 


.00000000 

-.00540841 

.02087231 

.05153854 

.24624152 

.23145357 

.03746083 

.00850816 


.00000000 

.00360077 

.01335611 

.02840036 

.05570463 

.21808803 

.17627681 

.01628154 


.00000000 

-.00245244 

.00866234 

.01742762 

.02036150 

-.04873011 

.15151413 

.08078532 


.00000000 

.00004205 

-.00330210 

.00651400 

-.01050027 

.01555127 

-.02205750 

. 04047503 . 
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- 1.00000000 -. 80075^0 -. 67718^8 
.00000000 .00000000 .00000000 
.03858768 .07025101 —.01272820 

-.01821052 .11002103 .14050210 

.01277058 -.03008507 .17700013 

-.00015870 .02555263 -.04830068 

.00630252 -.01682355 .02736602 

-.00303003 .01027426 -.01570330 

.00100655 -.00516256 .00765877 

-.00057030 .00146755 -.00214812 


.36311746 .00000000 .36311746 

.00000000 .00000000 .00000000 

.00666755 —.00427369 .00268860 

.02983602 .01657357 -.01065025 

.16134243 —.04256122 .02366286 

.21134521 .21602068 —.04608886 

.04608886 .21602068 .21134521 

.02366286 —.04256122 .16134243 

.01065025 .01657357 -.02683602 

.00268860 -.00427360 .00666755 


.67718628 .80675800 1.00000000 

.00000000 .00000000 .00000000 

.00214812 .00146755 —.00057039 

.00765877 -.00516256 .00166655 

.01570330 .01027426 -.00363063 

.02736692 -.01682355 .00630252 

.04836668 .02555263 -.00615870 

.17760013 -.03608567 .01277058 

.14050210 .11002163 -.01821052 

.01272820 .07025101 .03858768 


(B.7) 


-1.00000000 

-.01053301 

-.73877386 

-.47702405 

—.16527806 

.16527806 

.47702405 

.73877386 

.01053301 

i.oooooooon 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.00000000 

.03003030 

.05643504 

-.01021815 

.00535041 

—.00344025 

.00243838 

-.00180282 

.00133052 

-.00003228 

.00036504 

-.01478117 

.00711306 

.11445424 

-.02425882 

.01351771 

—.00001030 

.00647053 

-.00474040 

.00327262 

-.00127744 

.01060201 

-.03317004 

.14707011 

.15061048 

-.03544673 

.02011031 

-.01340704 

.00054557 

-.00647836 

.00251255 

-.00700220 

.02203034 

-.04165330 

.18070787 

.18731172 

-.04253068 

.02416450- 

-.01504323 

.01047601 

-.00401505 

.00578828 

-.01543233 

.02503502 

-.04471061 

.10460760 

.10460760 

-.04471061 

.02503502 

-.01543233 

.00578828 

-.00401505 

.01047601 

-.01504323 

.02416450 

-.04253068 

.18731172 

.18070787- 

-.04165330 

.02203034 

-.00700220 

.00251255 

-.00647836 

.00054557 

-.01340704 

.02011031 

—.03544673 

.15061048 

.14707011 

-.03317004 

.01060201 

-.00127744 

.00327262 

-.00474040 

.00647053 

-.00001030 

.01351771 

-.02425882 

.11445424 

.00711306 

-.01478117 

.00036504 

-.00003228 

.00133052 

-.00180282 

.00243838 

-.00344025 

.00535041- 

-.01021815 

.05643504 

.03003030. 


(B.8) 


- 148 - 



Newton Integrating Matrices 


The following weighting matrices are repeated from Ref. [27]. For the Newton 
weighting matrices, the subscript n denotes the degree of the assumed polynomial 
upon which the integrating matrix is based. The parameter h is the step size. 



0 0 
1 0 
1 1 


O' 

0 

0 


LO . . . 0 1 iJ 


(B.9) 


-0 

5 

0 



0 

0 

.0 


0 0 0 . . . 0 * 

8-1 0 ... 0 

5 8-1 ... 0 

5 8-10 

0 5 8 -1 

0-18 5. 


(RIO) 



■ 0 

0 

0 

0 

0 




0- 


9 

19 

-5 

1 

0 



. . . 

0 


-1 

13 

13 

-1 

0 



. . . 

0 

h 

0 

-1 

13 

13 

-1 




0 

11 

; 








* 


0 

. . . 



-1 

13 

13 

-1 

0 


0 

... 



0 

-1 

13 

13 

-1 


. 0 

... 



0 

1 

-5 

19 

9. 


149 - 



^5 


- 0 

0 

0 

0 

0 

0 

0 - 

251 

646 

-264 

106 

-19 

0 

0 

-19 

346 

456 

-74 

11 

0 

0 

0 

-19 

346 

456 

-74 

11 

0 



0 

-19 

346 

456 

-74 

11 

0 

0 

0 

-19 

346 

456 

-74 

11 

0 ... 

0 

11 

-74 

456 

346 

-19 

0 ... 

0 

-19 

106 

-264 

646 

251 . 
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0 

0 

0 

0 

0 

0 

0 




0 - 

475 

1427 

-798 

482 

-173 

27 

0 



. . . 

0 

-27 

637 

1022 

-258 

77 

-11 

0 



. . . 

0 

11 

-93 

802 

802 

-93 

11 

0 



. . . 

0 

0 

11 

-93 

802 

802 

-93 

11 




0 

0 

... 



11 

-93 

802 

802 

-93 

11 

0 

0 

. . . 



0 

11 

-93 

802 

802 

-93 

11 

0 




0 

-11 

77 

-258 

1022 

637 

-27 

0 




0 

27 

-173 

482 

-798 

1427 

475 . 
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0 

0 

0 

0 

0 

0 

0 

0 

0 - 

10087 

65112 

—46481 

37504 

—20211 

6312 

—863 

0 

0 

—883 

25128 

46080 

— 18256 

7200 

—2088 

271 

0 

0 

271 

—2760 

30810 

37504 

— 8771 

1608 

—101 

0 

0 

0 

271 

—2760 

30810 

37504 

—6771 

1608 

—101 

0 




h 

60480 


0 


0 

0 


0 

0 


271 —2760 30818 37504 —8771 1608 —101 0 

0 271 —2760 30810 37504 —6771 1608 —101 

0 —101 1608 —6771 37504 30810 —2760 271 

0 271 —2088 7200 —16256 46080 25128 —863 

0 —863 6312 —20211 37504 —46461 65112 10087 . 
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W7 = 



0 

0 

0 

0 

0 

0 

0 

0 

0 


On 


36700 

130840 

— 121707 

123133 

—88547 

41400 

— 11351 

1375 

0 


0 


— 1375 

47700 

101340 

— 44707 

26883 

—11547 

2000 

—351 

0 


0 


351 

— 4183 

57627 

81603 

—20227 

7227 

— 1710 

101 

0 


0 


— 101 

1870 

—0531 

68323 

68323 

—0531 

1870 

— 101 

0 


0 

h 

0 

— 101 

1870 

—0531 

68323 

68323 

—0531 

1870 

— 101 


0 

120060 

0 


—101 

1870 

—0531 

68323 

68323 

—0531 

1870 

— 101 

0 


0 


0 

— 101 

1870 

—0531 

68323 

68323 

—0531 

1870 

— 101 


0 


0 

101 

— 1710 

7227 

—20227 

81603 

57627 

—4183 

351 


0 


0 

—351 

2000 

—11547 

26883 

— 44707 

101340 

47700 

—1375 


0 


0 

1375 

— 11351 

41400 

—88547 

123133 

— 121707 

130840 

36700 - 
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A Sample Gauss-Legendre Integrating Matrix 

The Gauss-Legendre integrating matrix presented here is an example of an 
integrating matrix that does not use end point nodes. Obviously, such matrices 
cannot be used for intervals where boundary conditions must be applied at the end 
points; instead, their intended use is for interior integration regions. As pointed out 
in Section 3.2, the merging process, which forms global weighting matrices from the 
summation of local weighting matrices, allows one to combine weighting matrices 
with and without endpoint nodes. 

Gauss-Legendre integrating matrices, which are based upon Legendre polyno- 
mials, are related to Gauss-Legendre quadrature. These matrices are also a close 
relative of the Jacobi integrating matrices since Legendre polynomials prove to be 
a special case of the more general Jacobi polynomials with a = ^ = 0. In fact, the 
general calculation procedure in Appendix A for Jacobi weighting matrices, when 
adapted to Legendre polynomials, can be used for the Gauss-Legendre weighting 
matrices. Note, however, that the Gauss-Legendre matrices, because of the lack 
of end points, are rectangular rather than square. To make the Gauss-Legendre 
matrices “conform” to other matrices with end points, the first and last columns of 
the matrix are “padded” with zeroes as shown in the example below; this padding 
allows the definition of consistent merging rules as noted in Section 3.2. 

The weighting matrix below, with two internal grid points, is for the normalized 
interval [—1,1]; it can be transformed to [0,1] by multiplying by the factor 1/2. 
Similar to the Jacobi matrices, the coordinate values x, which denote grid points 
on [—1, 1], are listed above the horizontal bar in the matrix. The grid points t on 
[0,1] can be obtained with the linear transformation t == 0.5(1 + x). The subscript 
n refers to the number of discretization intervals, with the integration being exact 
for all polynomials of degree <2n — 3. 




-.57735027 .57735027 


.00000000 

.00000000 

.00000000 

-.00000000 


.00000000 

-.07735027 

.57735027 

.50000000 


.00000000 

.50000000 

.57735027 

-.07735027 


.00000000 

.00000000 

.00000000 

. 00000000 - 


(R2) 
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Appendix C 


Composite Laminate Constitutive Equations 


The behavior of multi-ply laminated composites and composite skin box-beams 
can be predicted by developing relationships in the laminate axes between the 
applied loads and the resulting deflections. To develop such relationships, one begins 
with the stress-strain behavior of a single lamina of an orthotropic material. The 
lamina behavior is subsequently transformed from principal axes of the lamina to 
the reference axes for the multi-ply laminate. By then applying the assumptions of 
thin plate theory, expressions are developed from the properties of the constituent 
laminae relating force and moment resultants acting on the plate to midplane strains 
and plate curvatures. These expressions are the laminate constitutive equations. 
The following is intended only as a brief description of the process for obtaining 
the constitutive equations. More detailed presentations of these developments are 
given by Jones [25] and Ashton, Halpin, and Petit [56]. 

Assuming a state of plane stress, the simplified stress-strain relationships for a 
lamina of orthotropic material can be written as 


I”'! 


Qll Qi2 0 


fl ■ 


r 

. = 

Qi 2 Q 22 0 



(C.l) 

iri2; 

k 

. 0 0 ^ 66 . 

k 

712 

ib 


where 


158 - 



Qii 


El 

1 - 1 ' 121^21 


Qi 2 


V 12 E 2 

1 - V12V21 


^ 2 \E\ 

1 - I'i2l'21 


<922 = 


E2 

1 - V 12 V 21 


<966 = <^12» 


((7.2) 


and the subscript k denotes the kth layer in the laminate. The terms Q{j, which 
are referred to as reduced stiffnesses, are defined in the lamina principle axes (see 
Fig. C-1). 

If 6 is defined as the angle between the i-axis and the 1-axis (see Fig. C-2), 
then a transformation of both stress and strain components in Eq. (C.l) leads to 
the following stress-strain equations for the lamina in the laminate axis system: 




<9ii <9 i2 <9i6 



= 

^12 ^22 §26 




.<9 16 <926 <966. 

k 

. Ixy ) 


(C.3) 


where the Q^j are the transformed reduced stiffnesses, which have the definitions 
Qll = <9il cos^ 0 -I- 2(<5i 2 + 2(566)sin^ Ocos^ 6 + Q 22 sin^ 6 

^12 =(Qll +<922 “4(566)sin^^cos^tf+Ql2(sin^^ + cos^^) 

Q 22 = <9 1 1 ^ + 2(<9i 2 + 2Q66) sin^ ^ cos^ ^ H- <922 <^os^ 0 

{CA) 

<9i 6 = (<9ll - <9 i 2 - 2(966) sin ^ cos^ ^ + (Qi 2 - Q 22 + 2^66) sin^ ^ cos ^ 

Q 26 = (<9ll - <9 i 2 - 2<566)sin^^cos^ + (<5 i 2 - <922 + 2(5e6) sin ^ cos^ ^ 

§66 = (<9ll + <922 - 2<9 i 2 - 2<366)sin^^cos^^ + <566(sin^^ + cos'll?) . 


A convenient form exists for calculating the transformed reduced stiffnesses in 
terms of stiffness invariants. This invariant form, which was first developed by Tsai 
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and Pagano [57], is written as 


Qli — U 1 + U 2 cos2<? + Uz cos4^ 
Qi2 = + cos4^ 

Q 22 = U 1 —U 2 cos 26 + Uz cos i$ 

QlZ — ^t /2 sin 26 + Uz sin 46 
Q 26 ~ ^U 2 sin 26 — Uz sin 4^ 

A 

Qee = ^5 - Uz cos 46 


where the stiffness invariants are 

Ui = -(3Qii + 3 Q 22 + 2Qi2 + 4(566) 

U2 = 2^Qii ~ Q22) 

= g(<5ll + Q 22 — 2(5i2 — 4(566) (C'-6) 

U 4 = ^(Qll + <522 +6(5 i2 “ 4(566) 

U 5 = 3(^11 + <922 — 2 Qi2 + 4(566) • 


With the foregoing knowledge of the behavior of a single layer, classical lamina- 
tion theory for thin plates can now be used to determine the behavior of the 
laminate, or laminated box-beam. Classical lamination theory embodies a collec- 
tion of stress and deformation assumptions which constitute the familiar Kirchhofif 
hypothesis for plates. Under these assumptions, the strains for any point in the 
laminate can be written in terms of geometrical midplane displacements «o, vq, and 


- 156 - 



ivq as 



duo 

d^wo 


Cx — 

dx 


f 


dvQ 

d^WQ 


= 

dy 

5y2 

f 


duo 

5t;o 

d^-wo 

'yxjf — 

dy 


£.Z 

dxdy 


(C.7) 


Furthermore, the laminate strains given by Eq. (C.7) are equivalently expressed in 
terms of the midplane strains and plate curvatures as 


' ' 



/ Kx ' 


. = . 

4 


.Ixy. 





(C.8) 


If 'we now define the vectors 


O’ = {CTX 

•^y 

^xy}^ , 


= {e° 

4 


(C.9) 


Ky 

Kary} j 



then, with the aid of Eq. (C.8), the A:th layer stress-strain relationships in Eq. (C.3) 
can be expressed in the form 




(c-.io) 


Force and moment resultants offer a convenient way of expressing the relation- 
ships between internal and external loads of a laminated plate. With this in mind, 
it is convenient to make use of a force resultant vector n and a moment resultant 
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vector m which appear as 


n = {N^ Ny N^y}'^, fn={M» My . ((7.11) 


By assuming that the plate (or box-beam) has total thickness t and a number of 
layers nl, the force and moment resultants (i.e., force and moment per unit width) 
are defined by integrating the stresses in each lamina over the thickness, that is, 



((7.12) 


By substituting Eq. (C.IO) into Eq. (C.12) and noting that [^, £®, and k are not 
functions of z, one obtains the force and moment resultant expressions 

f f^k r^k \ 

n=^[^jfc<£®/ dz + K zdz\ ((7.13) 

Jb=i I -'^k-l J^k-\ ] 

and 

{ f^k f^k ) 

m=^[^jfc<£®/ z dz + K I z^ dz\ . (C'-14) 

k=i [ ''^k-l ■'^Jfc-1 } 

After carrying out the integrations, Eqs. (C.13-C.14) can be written as 


' ■ 


•All 

Ai2 

Ai6 

Bn 

B\2 

Bn' 

•eO 

Ny 


Ai2 

A 22 

A 26 

Bn 

B22 

B26 

4 

Nxy 


Ai6 

A 26 

Ab6 

Bn 

^26 

^66 

^xy 

Mx 


Bn 

Bi2 

^16 

Dn 

D\2 

Dn 

Kx 

My 


Bi2 

B22 

^26 

D\2 

D22 

D26 

Ky 

Mxy. 


■^16 

B26 

^66 

Dn 

D26 

^66J 

.^xy. 
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where 


NL 


^ ] {Qii)k i^k ^k—l) 

it==l 


^ NL 

^ k=\ 

^ NL 

Dij = — y~^(Qtv)fe ~ i) 
A:==l 


(C.16) 


In Eqs. (C.15) and (C.16), the Aij are extensional stiffnesses, the B{j are coupling 
stiffnesses, and the D,y are bending stiffnesses. For those laminates with layers 
arranged symmetrically about the geometric midplane of the plate, the terms 
will be zero. 

The integrating matrix approach requires knowledge of the compliance terms 
rather than the stiffness terms. The compliance terms are obtained by simply 
inverting the (6 X 6) matrix of stiffness terms appearing in Eq. (C.15). In symbolic 
form, the composite compliance relationship can be written as 



((7.17) 


In Eq. (C.17), the subscripts of the terms in the (3X3) submatrices will be designated 
t>y i,j = 1,2,3, rather than i,j — 1,2,6. 
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Appendix D 


Modified Strip Theory Aerodynamics 


The following nondimensional unsteady aerodynamic loads have evolved from 
the modified strip theory developed by Yates [58,59]. These Laplace transformed 
aerodynamic loads, which are valid for arbitrary motion, are made applicable to 
the Laplace domain by inclusion of either a generalized Theodorsen’s function or 
an appropriate rational approximation as discussed by Edwards [60-62]. The loads 
presented here are intended for use with the Laplace transformed fiutter equations 
(cf., Eqs. (2.21) and (4.30) ), and correspond to those load terms required by Eq. 
(2.28). The nondimensionalization of the aerodynamic loads is the same as that 
indicated for the p’s in Eq. (2.31). 

As in the theory presented by Yates, these aerodynamics employ a variable sec- 
tion lift-curve slope ao(x) instead of 27 t, a variable section aerodynamic center adx) 
instead of —0.5 (the quarter-chord), and a variable section structural reference axis 
location a{x). The Theodorsen function can be modified by a factor which accounts 
for compressibility effects on the magnitude of the lift and pitching moment. For 
more details on the modifications accompanying a particular planform and Mach 
number, see the above references by Yates. It should be noted that for incompres- 
sible flow over untapered wings having a lift-curve slope of 2 jt and an aerodynamic 
center at the quarter-chord, the aerodynamics given here reduce to the elementary 
strip theory used by Barmby, Cunningham, and Garrick [63]. 

The nondimensional Laplace transformed aerodynamic loads providing lift and 
pitching moment per unit span can be conveniently calculated with the aid of an 
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aerodynamic influence matrix Qfc(«*, X). The terms appearing in Qjcd are given by 


Ia(a*,X) 
L'j{s*, X) 

Lr{s*,\) 

M4s*,X) 

Ma(«*,X) 

From the 
the form 


= — — RC(s)ao\/X 

= — 7T\/X ^s* — RC(5)ao(a(/ — a)\/X ^s* — RC(a)aox| 

= f|;r\/X /?s* + RC(s)aox| tan A 
= p s* — RC(s)ao(a<f — a)x| tan A 

= 6|;ray5^s*^ + RC(5)oo(a — ac)\/X/3s*| 

= 6^^{— n-(g + — ;r(a^ — a)\/\0s* + RC(6)oo(a — ~ a)\/\^s* 

+ RC(s)ao(a — ac)x| 

= 6^^|— 7ra\/X — RC(s)ao(a — ac)x| tan A 
= |~^(8 0s* + RC(5)oo(o — ac){aj — o)x|tanA. 

(P.l) 

terms in Eq. (D.l), the aerodynamic matrix Qf£>(«*>X) is constructed in 


Qf£>(«*)X) 


0 0 0 0 

0 0 0 0 

0 i.'f L.U) ^a,T 

0 'My Ma,r. 


{D.2) 


where 


and 


J^a,r ='Ma+'>^r-D, (DA) 

with D being a differentiating matrix. (Note that the subscript t is used as a 
reference to the variable, r == ||.) For moderate angles of wing aerodynamic sweep 
and large aspect ratio, the last term on the right-hand side of Eqs. (D.3-D.4) is 
negligible compared to the remaining term. 

The following definitions apply to the foregoing equations, in which X is a 
dimensionless dynamic pressure parameter and s* is a dimensionless Laplace vari- 
able: 


8 


* 



sb ^8* 


(EI)r 
y* = ycosA 


rriR 





(D.b) 


The aerodynamic downwash point is calculated from 


ao . 


(C.6) 


A convenient rational approximation to the Theodorsen function, attributed to 
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R.T. Jones, is given by 

, 0 . 552 + 0.28085 + 0.01365 

C(5) = — . {D.7) 

52 + 0.34555 + 0.01365 

As discussed in Yates [58], for higher Mach numbers M, C(5) can be empirically 
corrected for compressibility effects by scaling by a factor R(5, M), which is deter- 
mined from a ratio of the magnitudes of compressible and incompressible circulation 
functions. 
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Appendix E 


Solutions for Problems with 
Concentrated Loads 


Concentrated loads present a special type of discontinuity that can be treated 
with the help of delta functions. To include the possibility of delta functions in 
the integrating matrix approach, it is first necessary to extend the developments 
of Section 3.1 for the integrating matrices of continuous integrands. Once delta 
functions have been included in the formulation of integrating matrices, the solution 
methods presented in Chapter 4 can be expanded to accommodate concentrated 
loads. 

Instead of considering a simple continuous function f{x), as was done in Section 
3.1, it is now necessary to consider a function of the form 

7i^) = f{^) + Pi (^ 1 ) 

where pf gives the magnitude of the delta function (concentrated load) at the point 
X = Xi. Taking the integral of Eq. (E.l) over the subinterval [x^, x^^i], with the delta 
function located at yields 


r^i+i /•^t+i ^ y^t+i 

/ J{x)dx-= f{x)dx + pf^^ 6{x-xi+i)dx 

Jx- Jxi 

= / f{x)dx + pf^i. 

J X2 


(E.2) 


By using the definitions for {7} and {/} given in Eqs. (3.5) and (3.6) and including 
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the definition of the column vector 


} = (po » > Pjv)^» (£^-3) 

one can express the set of all integrals given by Eq. (E.2) as 

{n = '^n{f} + {p^}- (EA) 

This result is analogous to Eq. (3.7); note, however, that the first term in {T) is 
zero only if is zero. As in Section 3.1, Eq. (E.4) can be premultiplied by the 
summing matrix $ to yield the integrating matrix relationship 

(7}-M/} + S{p+} (E.h) 

where the integrating matrix L is the same as defined in Eq. (3.11). Thus, from 
Eq. (E.5) it is clear that the “integrating matrix” for delta functions is simply the 
summing matrix S. 

The foregoing results can be applied to solving problems as discussed in Section 
4.1. Consider, for example, the discretized equations given in Eq. (4.2). After 
including the concentrated loads, the equation can be written in the form 

y' = Zy - X(A + K'^6{x - *t))y - (sr + a;J‘^(a: - Xj)) . {E.d) 

Applying the integrating matrices then yields 

y = L(Zy - XAy - ir) - 5(XA+y + a+) + k. {E.7) 

This result compares with that of Eq. (4.3). Note that the global summing matrix S 
will be composed of diagonal blocks of dimension (AT+l)X(iV+l), which corresponds 



If one now solves for k, the equivalent of Eq. (4.9) is obtained. Thus, 

k = -B[L(Zy - XAy - a^) - S(XA+y + i+)] - B„/.y. {E.8) 

After substituting Eq. (E.8) into Eq. (E.7), grouping similar terms, and then rear- 
ranging, one has the result 


[H - X(FA + F+A+)]y = f (E.9) 

where 

H = I + B„fc + FZ 

F = [B-I]L, F+ = [B-I]S 

f = Far + F+i+. 

It is clear that these equations are simply an extension of Eqs. (4.10-4.13). 

Carrying the analysis one step further, the reduced nonhomogeneous linear 
system in Eq. (4.26), for Hpo — 0) appears as 

[I - X(TAf^o + T+A+p)]y^ = + T+5+ (E.13) 


(£’. 10 ) 

(£. 11 ) 

(£. 12 ) 


where 


T = (£.14) 

and 

T+ = ■ (^15) 

Similarly, the reduced eigenvalue problem can be written as 

[(TAfc 4- T+A+o) - (l/X)I]yp = 0 . (£.16) 


- 167 - 




Appendix F 


Constraint Equations 


Auxiliary constraint equations can be used in specific instances to reduce the 
degrees of freedom required for a numerical solution. To carry out reductions on 
matrix equations, it is convenient to think in terms of partitioned matrices, where 
the variables to be eliininated are placed in a partition of the solution vector, and 
all matrices appearing in the matrix equation are partitioned accordingly. The 
constraint equations are then used to obtain a transformation of variables 


y = Uy* 


where y is a vector containing the original variables, y* is a vector of transformed 
variables, and U is a transformation matrix. It is assumed initially that U is a 
square matrix that possesses an inverse, but as noted below, one can make use of 
other types of transformations that are characterized by rectangular transformation 
matrices. 

Consider first a partitioned system of linear equations 


Gil Gi2ifyi 


G21 G22, 


iD-i::i 


Equation (F.l) can be substituted into Eq. (F.2), and the resulting matrix equation 
premultiplied by If U does in fact represent an elimination of variables, then G 
will be transformed by this similarity transformation into a block triangular form 
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I 


such that 


■G'll 

G21 



where contains the eliminated variables and f* = U ^f. The reduced linear 
system is obtained directly from Eq. (F.3) as 


G'liyJ = fl 


{FA) 


An identical approach can be followed for the eigenvalue problem given by 


Gll 

Gi 2 jy 

'l-xD 

.G21 

G 22 Jly 

2/ ly2j 


The similarity transformation in this situation yields the equations 



from which the reduced eigenvalue problem is 


(F.5) 


(F.6) 


[Gll-^J]yi =0. (F.7) 

For those problems for which U is a square, invertable matrix, the foregoing 
reduction process is carried out as presented. On the other hand, it is also possible 
to use a transformation that involves a rectangular transformation matrix. In this 
situation, the above procedure requires a congruence transformation rather than a 
similarity transformation. To convert the resulting eigenvalue problem to standard 
form requires the inversion of the matrix product U^U. Unless this product results in 
a diagonal matrix, or otherwise has some special form, the best numerical approach 
makes use of singular value decomposition and the pseudoinverse. Discussions on 
this subject can be found in Strang [64] and Atkinson [65]. 
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Appendix G 


Calculation of Transition and 
Influence Matrices 


Integrating matrices ofifer a convenient method for numerically calculating tran- 
sition matrix solutions of two-point boundary problems. For structural solutions, 
the transition matrix can also be applied to the development of stiffness influence 
matrices. In formulating eigenvalue problems, the direct integrating matrix ap- 
proach presented in Chapter 4 is to be preferred for its simplicity and numeri- 
cal efficiency. Nevertheless, it is possible to use the transition matrix concept to 
advantage in certain types of numerical calculations. A detailed description of 
transition (also known as transfer or transmission) matrix methods in structural 
mechanics is presented by Pestel and Leckie [66]. A brief, but useful account of the 
method can also be found in Chapter 10 of McGuire and Gallagher [20]. In addition, 
the reader can refer to Chapter 7 of Boyce and DiPrima [67] for a general review 
of fundamental matrices and the role they play in the solution of linear differential 
equations. For useful applications and properties of transition matrices, one should 
consult Appendix A4 of Bryson and Ho [21] and Chapter 9 of Kailath [68]. 

The calculation of a transition matrix begins with the integrated version of the 
homogeneous state vector equations. These equations can be obtained from Eq. 
(4.3) by dropping the nonhomogeneous term Lir, thus yielding 

y = LZy-XLAy-l-k. (G.l) 

It will be assumed that the ordering of the components of the global state vector is 
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the same as that used in the derivations presented in Section 4.1. As will be recalled, 
the global state vector is arranged such that it can be partitioned into generalized 
force and generalized displacement subsets, and the discrete set of values for a 
particular variable are grouped together. (An alternate ordering scheme that can be 
used for transition matrix derivation involves consecutively placing each local state 
vector into the global state vector. This ordering scheme, however, requires one to 
use a modified definition of the global integrating matrix. The modified integrating 
matrix is obtained by expanding each element of L into a diagonal submatrix, with 
the element value repeated in each of the diagonal terms of the submatrix.) 

As the next step in calculating the transition matrix, choose the constant of 
integration to be equal to the local state vector at the end point a: = 0. This means 
that the global constant vector of integration can be written as 

'Bff Bfd~\ fy ^0 \ 

k = syo = { } (G.2) 

.Bdf Sddj ly ^0 ) 

where B will be termed a selection matrix since it selects the component of the local 
state vector to be placed in each component of k. B consists primarily of zeroes, 
but has appropriately placed unit terms. For clarity, it is worth noting that if NS 
is the number of state variables and {N + 1) is the number of grid points, then the 
dimension of yo will be NS and the dimensions of B will be NS(AT + 1) XNS. 

If one now combines Eq. (G.l) with Eq. (G.2) the result is 

Hy = Byo = k (G.3) 

where 

H = 1- 

The corresponding partitioned form of 

.H/jf' Wdd\ 




In a straight forward manner, Eq. (G.5) can be solved by partitions, which yields 


where 


Tff 

Yfd 

"^DF 

Td/) 


fy^l [Tff TF/jiryFoI 

lyp/ LTdf T£)£>Jlyr»oJ 


[Sff — HfdTdf] 

HJp[Bfd — HFDTx>i>j 

[Hz>£» — H£»FH7/HFi)]“^[S£)F — HdfH^jJHff] 
(H/jz) — H/jfHff^f/?]~^(Sc>z? — Hdf^ff^fd] ■ 


(G.6) 


(G.7) 


By definition, the transition (transfer) matrix transforms the state vector at 
one point into the state vector at another point. Considering the transition matrix 
between the state vector at x = 0 and the state vector at the ith grid point, one 
can write the general transition matrix relationship as 

■ ^FF{ ^FDi I f ' 

At the same time, it is possible to write another expression for the local state vector 
appearing on the left-hand side of Eq. (G.8). This expression, which makes use of 
another selection matrix, F,-, transforms the global state vector into the local vector 
at the »th point. The expression is written as 


fy^i I _ [ ^FDi j (YFo \ 

\yi>J ^£T>Jly/>oi 


(G.8) 


D-i:' ;i;i 


(G.9) 


By substituting the result for the global state vector from Eq. (G.6) into Eq. (G.9), 
and then comparing with Eq. (G.8), one finds that the transition matrix for the tth 
grid point is given by 

TjTFF 

riToF riTDD. ' 




(G.IO) 
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This equation, in conjunction with Eq. (G.7), provides a method for calculating the 
transition matrix. 

The foregoing results for transition matrices can be applied in calculating a 
stiffness influence matrix. For the influence matrix calculation, the load terms 
multiplied by X in Eq. (G.l) are set to zero; for certain problems this can mean 
that Hfd in Eqs. (G.5) and (G.7) will be zero. The strategy in developing a 
stiffness influence matrix from a transition matrix is to obtain a force-displacement 
relationship between degrees of freedom at each end of the normalized interval 
[0,1]. By definition, the matrix that relates the force degrees of freedom to the 
displacement degrees of freedom will be the influence matrix. 

By making use of Eq. (G.8), and taking into account the properties of a 
normalized transition matrix, one can write the expression 

fy/>(0)| r^2>F(0) ^/>rj(0)-||yFo| 

lyr»(l)i L^z)f(I) ^o£)(l)Jly^o/ 

Solving for the constant vector on the far right-hand side of Eq. (G.ll) gives 


Similar to Eq. (G.ll), an expression can be written for the forces at the 
endpoints. This expression appears as 

JyF(o)i r4^FF(o) 4>FZj(o)ijyFoi r i o -|fyFoi 

l~yF(l)/ [— ^Ff(I) — 4^Fd( 1)J ly^o/ L— ^Ff(I) — ^Fd( 1)J ly^o/ 

Substituting the constant vector from Eq. (G.12) into Eq. (G.13) yields 



1 0 


|yr>(0)] 

yz?(i)J 


Jl 


(C?.12) 


.^of(1) <Pdo(1), 


I 

o(l)JW/ 


(G.ll) 



where the stiffness influence matrix K' is given by 


K' 




(G.15) 


Provided that certain conditions are met by the transition matrix, it can be 
shown that K' is a symmetric matrix. To demonstrate this fact requires knowledge of 
the transition matrix for the adjoint system of homogeneous differential equations, 
where the adjoint equations are obtained by substituting — for Z in Eq. (G.l). 
The transition matrices of the original and adjoint systems share the identity 


= I. 


(G.16) 


where is the transition matrix for the adjoint system. Normally, one would 
obtain by solving the adjoint differential equations, but for symplectic systems 
(see Section 2.1) the adjoint transition matrix can be obtained directly from the 
original transition matrix through the relationship 

(G.17) 

where J is as defined in Chapter 2 (see Bryson and Ho [21], p.l57). With the aid 
of this relationship, it can be shown that if the numerically calculated transition 
matrices satisfy Eq. (G.16), then K' — = 0, which shows that K' is symmetric. 
It should be noted, however, that sample numerical calculations seem to indicate 
that the discretization level determines how accurately the numerically calculated 
transition matrices match the identity in Eq. (G.16). The satisfaction of Eq. (G.16) 
possibly might serve as an indicator of sufficient discretization, but at the present 
time this has not been verified and thus remains an object for further study. 
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A final item to be noted about transition matrix solutions is that calculations 
for nonhomogeneous linear problems can be simplified considerably if they are 
symplectic. That is, by making use of integrating matrices and by applying the 
relationship in Eq. (G.17), one finds that the usual variation of parameters solution 
for the state vector equations can be written in a conveniently calculable form. 
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